From: bangerth Date: Fri, 25 Aug 2006 17:23:48 +0000 (+0000) Subject: Factor out common code. X-Git-Url: https://gitweb.dealii.org/cgi-bin/gitweb.cgi?a=commitdiff_plain;h=c788175b3f74b7590b0fc8a66eef845bf9e8012d;p=dealii-svn.git Factor out common code. git-svn-id: https://svn.dealii.org/trunk@13748 0785d39b-7218-0410-832d-ea1e28bc413d --- diff --git a/tests/fe/abf_projection_01.cc b/tests/fe/abf_projection_01.cc index d48ab21597..3b4decec2f 100644 --- a/tests/fe/abf_projection_01.cc +++ b/tests/fe/abf_projection_01.cc @@ -1,5 +1,5 @@ -//---------------------------- abf_projection_01.cc --------------------------- -// abf_projection_01.cc,v 1.3 2003/06/09 16:00:38 wolf Exp +//---------------------------- abfprojection_01.cc --------------------------- +// abfprojection_01.cc,v 1.3 2003/06/09 16:00:38 wolf Exp // Version: // // Copyright (C) 2003, 2005, 2006 by the deal.II authors @@ -9,7 +9,7 @@ // to the file deal.II/doc/license.html for the text and // further information on this license. // -//---------------------------- abf_projection_01.cc --------------------------- +//---------------------------- abfprojection_01.cc --------------------------- /* * Project the function [1,1] onto a deformed grid and see whether the ABF @@ -18,739 +18,12 @@ -#include "../tests.h" -#include +char logname[] = "abf_projection_01/output"; +#include "deformed_projection.h" -#define PRECISION 2 -#include - -std::ofstream logfile ("abf_projection_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 - - -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 &, - const unsigned int ) const -{ - return (1.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)); - - // Parabolic inflow profile - for (unsigned int iCount = 0; iCount < this->n_components; iCount++) - return_value (iCount) = value (p, iCount); -} - - -/* - * Check the value of the derivative field. - */ - -void EvaluateDerivative (DoFHandler<2> *dof_handler, - Vector &solution) -{ - // This quadrature rule determines the points, where the - // derivative will be evaluated. - QGauss<2> quad (3); - FEValues<2> fe_values (dof_handler->get_fe (), quad, - UpdateFlags(update_values | - update_q_points | - update_gradients | - update_JxW_values)); - - const unsigned int n_q_points = quad.n_quadrature_points; - const unsigned int n_components = dof_handler->get_fe().n_components(); - const unsigned int dofs_per_cell = dof_handler->get_fe().dofs_per_cell; - - // Cell iterators - DoFHandler<2>::active_cell_iterator cell = dof_handler->begin_active(), - endc = dof_handler->end(); - - double err_l2 = 0, - err_hdiv = 0; - - std::vector local_dof_indices (dofs_per_cell); - - for (; cell!=endc; ++cell) - { - cell->get_dof_indices (local_dof_indices); - - fe_values.reinit (cell); - - // Get function values - std::vector > this_value(n_q_points, - Vector(n_components)); - fe_values.get_function_values (solution, this_value); - - // Get values from solution vector (For Trap.Rule) - std::vector > > - grads_here (n_q_points, std::vector > (n_components)); - fe_values.get_function_grads (solution, grads_here); - - for (unsigned int q_point=0; q_point -void create_mass_matrix (const Mapping &mapping, - const DoFHandler &dof, - const Quadrature &q, - SparseMatrix &matrix, - const Function &rhs_function, - Vector &rhs_vector, - const Function * const coefficient = 0) +void test () { - UpdateFlags update_flags = UpdateFlags(update_values | update_JxW_values | update_q_points); - if (coefficient != 0) - update_flags = UpdateFlags (update_flags | update_q_points); - - FEValues fe_values (mapping, dof.get_fe(), q, update_flags); - - const unsigned int dofs_per_cell = fe_values.dofs_per_cell, - n_q_points = fe_values.n_quadrature_points; - const FiniteElement &fe = fe_values.get_fe(); - const unsigned int n_components = fe.n_components(); - - Assert(coefficient == 0 || - coefficient->n_components==1 || - coefficient->n_components==n_components, ExcInternalError()); - - FullMatrix cell_matrix (dofs_per_cell, dofs_per_cell); - Vector cell_vector (dofs_per_cell); - std::vector coefficient_values (n_q_points); - std::vector > coefficient_vector_values (n_q_points, - Vector (n_components)); - - std::vector dof_indices (dofs_per_cell); - - std::vector > rhs_values(n_q_points, Vector(n_components)); - - typename DoFHandler::active_cell_iterator cell = dof.begin_active (), - endc = dof.end (); - for (; cell!=endc; ++cell) - { - fe_values.reinit (cell); - - cell_matrix = 0; - cell->get_dof_indices (dof_indices); - - const std::vector &weights = fe_values.get_JxW_values (); - rhs_function.vector_value_list (fe_values.get_quadrature_points(), rhs_values); - cell_vector = 0; - - if (coefficient != 0) - { - if (coefficient->n_components==1) - { - coefficient->value_list (fe_values.get_quadrature_points(), - coefficient_values); - for (unsigned int point=0; pointvector_value_list (fe_values.get_quadrature_points(), - coefficient_vector_values); - for (unsigned int point=0; point sign_change (dofs_per_cell, 1.0); - const unsigned int dofs_per_face = fe.dofs_per_face; - std::vector face_dof_indices (dofs_per_face); - - for (unsigned int f = 0; f < 2; ++f) - { - typename DoFHandler::active_face_iterator face = cell->face (f); - if (!face->at_boundary ()) - { - unsigned int nn = cell->neighbor_of_neighbor (f); - sprintf (buf, "Face %i NeigNeig %i\n", f, nn); - deallog << buf; - if (nn > 1) - { - face->get_dof_indices (face_dof_indices); - for (unsigned int j = 0; j < dofs_per_face; ++j) - { - sign_change[face_dof_indices[j]] = -1.0; - sprintf (buf, "DoF %i\n", face_dof_indices[j]); - deallog << buf; - } - } - } - } - - for (unsigned int point=0; point > val_vector (dofs_per_cell, - Vector (n_components)); - - // Precompute the component values - for (unsigned int i=0; i < dofs_per_cell; ++i) - for (unsigned int comp_i = 0; comp_i < fe.n_components (); - ++comp_i) - { - val_vector[i](comp_i) = sign_change[i] * - fe_values.shape_value_component(i,point,comp_i); - } - // Now eventually switch the sign of some of the ansatzfunctions. - // TODO - - for (unsigned int i=0; i -void create_right_hand_side (const Mapping &mapping, - const DoFHandler &dof_handler, - const Quadrature &quadrature, - const Function &rhs_function, - Vector &rhs_vector) -{ - const FiniteElement &fe = dof_handler.get_fe(); - Assert (fe.n_components() == rhs_function.n_components, - ExcInternalError()); - Assert (rhs_vector.size() == dof_handler.n_dofs(), - ExcDimensionMismatch(rhs_vector.size(), dof_handler.n_dofs())); - rhs_vector = 0; - - UpdateFlags update_flags = UpdateFlags(update_values | - update_q_points | - update_JxW_values); - FEValues fe_values (mapping, fe, quadrature, update_flags); - - const unsigned int dofs_per_cell = fe_values.dofs_per_cell, - n_q_points = fe_values.n_quadrature_points, - n_components = fe.n_components(); - - std::vector dofs (dofs_per_cell); - Vector cell_vector (dofs_per_cell); - - typename DoFHandler::active_cell_iterator cell = dof_handler.begin_active(), - endc = dof_handler.end(); - - if (n_components==1) - { - std::vector rhs_values(n_q_points); - - for (; cell!=endc; ++cell) - { - fe_values.reinit(cell); - - const std::vector &weights = fe_values.get_JxW_values (); - rhs_function.value_list (fe_values.get_quadrature_points(), rhs_values); - - cell_vector = 0; - for (unsigned int point=0; pointget_dof_indices (dofs); - - for (unsigned int i=0; i > rhs_values(n_q_points, Vector(n_components)); - - for (; cell!=endc; ++cell) - { - fe_values.reinit(cell); - - const std::vector &weights = fe_values.get_JxW_values (); - rhs_function.vector_value_list (fe_values.get_quadrature_points(), rhs_values); - - cell_vector = 0; - for (unsigned int point=0; pointget_dof_indices (dofs); - - for (unsigned int i=0; i -void project (const Mapping &mapping, - const DoFHandler &dof, - const ConstraintMatrix &constraints, - const Quadrature &quadrature, - const Function &function, - Vector &vec, - const bool enforce_zero_boundary = false, - const Quadrature & = QGauss2(), - const bool project_to_boundary_first = false) -{ - Assert (dof.get_fe().n_components() == function.n_components, - ExcInternalError()); - - const FiniteElement &fe = dof.get_fe(); - - // make up boundary values - std::map boundary_values; - - if (enforce_zero_boundary == true) - // no need to project boundary - // values, but enforce - // homogeneous boundary values - // anyway - { - // loop over all boundary faces - // to get all dof indices of - // dofs on the boundary. note - // that in 3d there are cases - // where a face is not at the - // boundary, yet one of its - // lines is, and we should - // consider the degrees of - // freedom on it as boundary - // nodes. likewise, in 2d and - // 3d there are cases where a - // cell is only at the boundary - // by one vertex. nevertheless, - // since we do not support - // boundaries with dimension - // less or equal to dim-2, each - // such boundary dof is also - // found from some other face - // that is actually wholly on - // the boundary, not only by - // one line or one vertex - typename DoFHandler::active_face_iterator face = dof.begin_active_face(), - endf = dof.end_face(); - std::vector face_dof_indices (fe.dofs_per_face); - for (; face!=endf; ++face) - if (face->at_boundary()) - { - face->get_dof_indices (face_dof_indices); - for (unsigned int i=0; i::type boundary_functions; - for (unsigned char c=0; c<255; ++c) - boundary_functions[c] = &function; - project_boundary_values (dof, boundary_functions, q_boundary, - boundary_values); -*/ - } - - - // set up mass matrix and right hand side - vec.reinit (dof.n_dofs()); - SparsityPattern sparsity(dof.n_dofs(), - dof.n_dofs(), - dof.max_couplings_between_dofs()); - DoFTools::make_sparsity_pattern (dof, sparsity); - constraints.condense (sparsity); - - SparseMatrix mass_matrix (sparsity); - Vector tmp (mass_matrix.n()); - - create_mass_matrix (mapping, dof, quadrature, mass_matrix, function, tmp); - sprintf (buf, "MM created\n"); - deallog << buf; - create_right_hand_side (mapping, dof, quadrature, function, tmp); - sprintf (buf, "RHS created\n"); - deallog << buf; - - constraints.condense (mass_matrix); - constraints.condense (tmp); - if (boundary_values.size() != 0) - MatrixTools::apply_boundary_values (boundary_values, - mass_matrix, vec, tmp, - true); - - SolverControl control(1000,1e-16); - PrimitiveVectorMemory<> memory; - SolverCG<> cg(control,memory); - - PreconditionSSOR<> prec; - prec.initialize(mass_matrix, 1.2); - // solve - cg.solve (mass_matrix, vec, tmp, prec); - - // distribute solution - constraints.distribute (vec); -} - - -int create_tria (unsigned int elm, Triangulation<2> &tria) -{ - std::vector > points_glob; - std::vector > points; - - points_glob.push_back (Point<2> (0.0, 0.0)); - points_glob.push_back (Point<2> (1.0, 0.0)); - points_glob.push_back (Point<2> (1.0, 0.5)); - points_glob.push_back (Point<2> (1.0, 1.0)); - points_glob.push_back (Point<2> (0.6, 0.5)); - points_glob.push_back (Point<2> (0.5, 1.0)); - points_glob.push_back (Point<2> (0.0, 1.0)); - - // Prepare cell data - std::vector > cells (1); - - switch (elm) - { - case 0: - cells[0].vertices[0] = 0; - cells[0].vertices[1] = 1; - cells[0].vertices[2] = 4; - cells[0].vertices[3] = 2; - cells[0].material_id = 0; - break; - case 1: - cells[0].vertices[0] = 4; - cells[0].vertices[1] = 2; - cells[0].vertices[2] = 5; - cells[0].vertices[3] = 3; - cells[0].material_id = 0; - break; - case 2: - cells[0].vertices[0] = 0; - cells[0].vertices[1] = 4; - cells[0].vertices[2] = 6; - cells[0].vertices[3] = 5; - cells[0].material_id = 0; - break; - } - - for (unsigned int i = 0; i < 4; ++i) - { - points.push_back (points_glob[cells[0].vertices[i]]); - cells[0].vertices[i] = i; - } - - tria.create_triangulation (points, cells, SubCellData()); - - return (0); -} - - -int plot_shapes (DoFHandler<2> &dof_handler) -{ - Vector solution (dof_handler.n_dofs ()); - std::set face_dofs; - - // Create set of all DoFs, which are on the boundary. - DoFHandler<2>::active_face_iterator face = dof_handler.begin_active_face(), - endf = dof_handler.end_face(); - std::vector face_dof_indices (dof_handler.get_fe().dofs_per_face); - for (; face!=endf; ++face) - { - face->get_dof_indices (face_dof_indices); - for (unsigned int i=0; i::iterator face_dof_iter; - for (face_dof_iter = face_dofs.begin (); - face_dof_iter != face_dofs.end (); ++face_dof_iter) - { - unsigned int dof = *face_dof_iter; - - sprintf (buf, "%i\n", dof); - deallog << buf; - - solution (dof) = 1.0; - - // Test the core functionality2 - 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 (4, 0); - - data_out->write_gnuplot (deallog.get_file_stream()); - - delete data_out; - - solution (dof) = 0.0; - } - - return (0); -} - - -int main (int /*argc*/, char **/*argv*/) -{ - 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; - FE_ABF<2> fe (0); - deallog << "Dofs/cell " << fe.dofs_per_cell - << "Dofs/face " << fe.dofs_per_face << std::endl; - - for (unsigned int elm = 0; elm < 3; ++elm) - { - create_tria (elm, tria_test); - - // Create a DoFHandler - dof_handler = new DoFHandler<2> (tria_test); - dof_handler->distribute_dofs (fe); - - deallog << "Dofs total " << dof_handler->n_dofs () << std::endl; - - // Alloc some DoFs - Vector solution; - solution.reinit (dof_handler->n_dofs ()); - - // Project solution onto FE field - ConstraintMatrix hn_constraints; - hn_constraints.clear (); - DoFTools::make_hanging_node_constraints (*dof_handler, - hn_constraints); - hn_constraints.close (); - - MappingQ1<2> map_default; - - project (map_default, *dof_handler, hn_constraints, - QGauss6<2> (), TestMap1<2>(2), - solution); - - // 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 (4); - - data_out->write_gnuplot (deallog.get_file_stream()); - - // Now write face DoFs .. - for (unsigned int i = 0; i < 4; ++i) - { - sprintf (buf, "%i %e\n", i, solution (i)); - deallog << buf; - } - - // Clean up ... - delete data_out; - delete (dof_handler); - tria_test.clear (); - } - - return (0); + check (fe); } diff --git a/tests/fe/abf_projection_01/cmp/generic b/tests/fe/abf_projection_01/cmp/generic index 84f06527b8..0f500d410f 100644 --- a/tests/fe/abf_projection_01/cmp/generic +++ b/tests/fe/abf_projection_01/cmp/generic @@ -2,8 +2,8 @@ DEAL::Dofs/cell 6Dofs/face 1 DEAL::Dofs total 6 DEAL::MM created -RHS created -Starting value 1.99 +DEAL::RHS created +DEAL:cg::Starting value 1.99 DEAL:cg::Convergence step 9 value 0 # This file was generated by the deal.II library. @@ -43,14 +43,14 @@ DEAL:cg::Convergence step 9 value 0 1.00 0.50 1.00 1.00 -DEAL::0 -1.000000e-01 -1 5.000000e-01 -2 1.000000e+00 -3 4.000000e-01 -Dofs total 6 +DEAL::-0.10 +DEAL::0.50 +DEAL::1.00 +DEAL::0.40 +DEAL::Dofs total 6 DEAL::MM created -RHS created -Starting value 1.19 +DEAL::RHS created +DEAL:cg::Starting value 1.19 DEAL:cg::Convergence step 8 value 0 # This file was generated by the deal.II library. @@ -90,14 +90,14 @@ DEAL:cg::Convergence step 8 value 0 1.00 1.00 1.00 1.00 -DEAL::0 6.000000e-01 -1 5.000000e-01 -2 4.000000e-01 -3 5.000000e-01 -Dofs total 6 +DEAL::0.60 +DEAL::0.50 +DEAL::0.40 +DEAL::0.50 +DEAL::Dofs total 6 DEAL::MM created -RHS created -Starting value 1.99 +DEAL::RHS created +DEAL:cg::Starting value 1.99 DEAL:cg::Convergence step 9 value 0 # This file was generated by the deal.II library. @@ -137,7 +137,7 @@ DEAL:cg::Convergence step 9 value 0 0.50 1.00 1.00 1.00 -DEAL::0 1.000000e+00 -1 6.000000e-01 -2 1.000000e-01 -3 5.000000e-01 +DEAL::1.00 +DEAL::0.60 +DEAL::0.10 +DEAL::0.50 diff --git a/tests/fe/deformed_projection.h b/tests/fe/deformed_projection.h new file mode 100644 index 0000000000..03e6274d71 --- /dev/null +++ b/tests/fe/deformed_projection.h @@ -0,0 +1,753 @@ +//---------------------------- deformed_projection.h --------------------------- +// deformed_projection.h,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. +// +//---------------------------- deformed_projection.h --------------------------- + +/* + * Project the function [1,1] onto a deformed grid and see whether the element + * elements can represent it exactly. This is the common code for the various + * *_projection_* tests. + */ + + + +#include "../tests.h" +#include + +#define PRECISION 2 + +#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 + +void test(); + + +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 &, + const unsigned int ) const +{ + return (1.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)); + + // Parabolic inflow profile + for (unsigned int iCount = 0; iCount < this->n_components; iCount++) + return_value (iCount) = value (p, iCount); +} + + +/* + * Check the value of the derivative field. + */ + +void EvaluateDerivative (DoFHandler<2> *dof_handler, + Vector &solution) +{ + // This quadrature rule determines the points, where the + // derivative will be evaluated. + QGauss<2> quad (3); + FEValues<2> fe_values (dof_handler->get_fe (), quad, + UpdateFlags(update_values | + update_q_points | + update_gradients | + update_JxW_values)); + + const unsigned int n_q_points = quad.n_quadrature_points; + const unsigned int n_components = dof_handler->get_fe().n_components(); + const unsigned int dofs_per_cell = dof_handler->get_fe().dofs_per_cell; + + // Cell iterators + DoFHandler<2>::active_cell_iterator cell = dof_handler->begin_active(), + endc = dof_handler->end(); + + double err_l2 = 0, + err_hdiv = 0; + + std::vector local_dof_indices (dofs_per_cell); + + for (; cell!=endc; ++cell) + { + cell->get_dof_indices (local_dof_indices); + + fe_values.reinit (cell); + + // Get function values + std::vector > this_value(n_q_points, + Vector(n_components)); + fe_values.get_function_values (solution, this_value); + + // Get values from solution vector (For Trap.Rule) + std::vector > > + grads_here (n_q_points, std::vector > (n_components)); + fe_values.get_function_grads (solution, grads_here); + + for (unsigned int q_point=0; q_point +void create_mass_matrix (const Mapping &mapping, + const DoFHandler &dof, + const Quadrature &q, + SparseMatrix &matrix, + const Function &rhs_function, + Vector &rhs_vector, + const Function * const coefficient = 0) +{ + UpdateFlags update_flags = UpdateFlags(update_values | update_JxW_values | update_q_points); + if (coefficient != 0) + update_flags = UpdateFlags (update_flags | update_q_points); + + FEValues fe_values (mapping, dof.get_fe(), q, update_flags); + + const unsigned int dofs_per_cell = fe_values.dofs_per_cell, + n_q_points = fe_values.n_quadrature_points; + const FiniteElement &fe = fe_values.get_fe(); + const unsigned int n_components = fe.n_components(); + + Assert(coefficient == 0 || + coefficient->n_components==1 || + coefficient->n_components==n_components, ExcInternalError()); + + FullMatrix cell_matrix (dofs_per_cell, dofs_per_cell); + Vector cell_vector (dofs_per_cell); + std::vector coefficient_values (n_q_points); + std::vector > coefficient_vector_values (n_q_points, + Vector (n_components)); + + std::vector dof_indices (dofs_per_cell); + + std::vector > rhs_values(n_q_points, Vector(n_components)); + + typename DoFHandler::active_cell_iterator cell = dof.begin_active (), + endc = dof.end (); + for (; cell!=endc; ++cell) + { + fe_values.reinit (cell); + + cell_matrix = 0; + cell->get_dof_indices (dof_indices); + + const std::vector &weights = fe_values.get_JxW_values (); + rhs_function.vector_value_list (fe_values.get_quadrature_points(), rhs_values); + cell_vector = 0; + + if (coefficient != 0) + { + if (coefficient->n_components==1) + { + coefficient->value_list (fe_values.get_quadrature_points(), + coefficient_values); + for (unsigned int point=0; pointvector_value_list (fe_values.get_quadrature_points(), + coefficient_vector_values); + for (unsigned int point=0; point sign_change (dofs_per_cell, 1.0); + const unsigned int dofs_per_face = fe.dofs_per_face; + std::vector face_dof_indices (dofs_per_face); + + for (unsigned int f = 0; f < 2; ++f) + { + typename DoFHandler::active_face_iterator face = cell->face (f); + if (!face->at_boundary ()) + { + unsigned int nn = cell->neighbor_of_neighbor (f); + deallog << "Face " << f << " NeigNeig " << nn << std::endl; + if (nn > 1) + { + face->get_dof_indices (face_dof_indices); + for (unsigned int j = 0; j < dofs_per_face; ++j) + { + sign_change[face_dof_indices[j]] = -1.0; + deallog << "DoF " << face_dof_indices[j] << std::endl; + } + } + } + } + + for (unsigned int point=0; point > val_vector (dofs_per_cell, + Vector (n_components)); + + // Precompute the component values + for (unsigned int i=0; i < dofs_per_cell; ++i) + for (unsigned int comp_i = 0; comp_i < fe.n_components (); + ++comp_i) + { + val_vector[i](comp_i) = sign_change[i] * + fe_values.shape_value_component(i,point,comp_i); + } + // Now eventually switch the sign of some of the ansatzfunctions. + // TODO + + for (unsigned int i=0; i +void create_right_hand_side (const Mapping &mapping, + const DoFHandler &dof_handler, + const Quadrature &quadrature, + const Function &rhs_function, + Vector &rhs_vector) +{ + const FiniteElement &fe = dof_handler.get_fe(); + Assert (fe.n_components() == rhs_function.n_components, + ExcInternalError()); + Assert (rhs_vector.size() == dof_handler.n_dofs(), + ExcDimensionMismatch(rhs_vector.size(), dof_handler.n_dofs())); + rhs_vector = 0; + + UpdateFlags update_flags = UpdateFlags(update_values | + update_q_points | + update_JxW_values); + FEValues fe_values (mapping, fe, quadrature, update_flags); + + const unsigned int dofs_per_cell = fe_values.dofs_per_cell, + n_q_points = fe_values.n_quadrature_points, + n_components = fe.n_components(); + + std::vector dofs (dofs_per_cell); + Vector cell_vector (dofs_per_cell); + + typename DoFHandler::active_cell_iterator cell = dof_handler.begin_active(), + endc = dof_handler.end(); + + if (n_components==1) + { + std::vector rhs_values(n_q_points); + + for (; cell!=endc; ++cell) + { + fe_values.reinit(cell); + + const std::vector &weights = fe_values.get_JxW_values (); + rhs_function.value_list (fe_values.get_quadrature_points(), rhs_values); + + cell_vector = 0; + for (unsigned int point=0; pointget_dof_indices (dofs); + + for (unsigned int i=0; i > rhs_values(n_q_points, Vector(n_components)); + + for (; cell!=endc; ++cell) + { + fe_values.reinit(cell); + + const std::vector &weights = fe_values.get_JxW_values (); + rhs_function.vector_value_list (fe_values.get_quadrature_points(), rhs_values); + + cell_vector = 0; + for (unsigned int point=0; pointget_dof_indices (dofs); + + for (unsigned int i=0; i +void project (const Mapping &mapping, + const DoFHandler &dof, + const ConstraintMatrix &constraints, + const Quadrature &quadrature, + const Function &function, + Vector &vec, + const bool enforce_zero_boundary = false, + const Quadrature & = QGauss2(), + const bool project_to_boundary_first = false) +{ + Assert (dof.get_fe().n_components() == function.n_components, + ExcInternalError()); + + const FiniteElement &fe = dof.get_fe(); + + // make up boundary values + std::map boundary_values; + + if (enforce_zero_boundary == true) + // no need to project boundary + // values, but enforce + // homogeneous boundary values + // anyway + { + // loop over all boundary faces + // to get all dof indices of + // dofs on the boundary. note + // that in 3d there are cases + // where a face is not at the + // boundary, yet one of its + // lines is, and we should + // consider the degrees of + // freedom on it as boundary + // nodes. likewise, in 2d and + // 3d there are cases where a + // cell is only at the boundary + // by one vertex. nevertheless, + // since we do not support + // boundaries with dimension + // less or equal to dim-2, each + // such boundary dof is also + // found from some other face + // that is actually wholly on + // the boundary, not only by + // one line or one vertex + typename DoFHandler::active_face_iterator face = dof.begin_active_face(), + endf = dof.end_face(); + std::vector face_dof_indices (fe.dofs_per_face); + for (; face!=endf; ++face) + if (face->at_boundary()) + { + face->get_dof_indices (face_dof_indices); + for (unsigned int i=0; i::type boundary_functions; + for (unsigned char c=0; c<255; ++c) + boundary_functions[c] = &function; + project_boundary_values (dof, boundary_functions, q_boundary, + boundary_values); +*/ + } + + + // set up mass matrix and right hand side + vec.reinit (dof.n_dofs()); + SparsityPattern sparsity(dof.n_dofs(), + dof.n_dofs(), + dof.max_couplings_between_dofs()); + DoFTools::make_sparsity_pattern (dof, sparsity); + constraints.condense (sparsity); + + SparseMatrix mass_matrix (sparsity); + Vector tmp (mass_matrix.n()); + + create_mass_matrix (mapping, dof, quadrature, mass_matrix, function, tmp); + deallog << "MM created" << std::endl; + create_right_hand_side (mapping, dof, quadrature, function, tmp); + deallog << "RHS created" << std::endl; + + constraints.condense (mass_matrix); + constraints.condense (tmp); + if (boundary_values.size() != 0) + MatrixTools::apply_boundary_values (boundary_values, + mass_matrix, vec, tmp, + true); + + SolverControl control(1000,1e-16); + PrimitiveVectorMemory<> memory; + SolverCG<> cg(control,memory); + + PreconditionSSOR<> prec; + prec.initialize(mass_matrix, 1.2); + // solve + cg.solve (mass_matrix, vec, tmp, prec); + + // distribute solution + constraints.distribute (vec); +} + + +int create_tria (unsigned int elm, Triangulation<2> &tria) +{ + std::vector > points_glob; + std::vector > points; + + points_glob.push_back (Point<2> (0.0, 0.0)); + points_glob.push_back (Point<2> (1.0, 0.0)); + points_glob.push_back (Point<2> (1.0, 0.5)); + points_glob.push_back (Point<2> (1.0, 1.0)); + points_glob.push_back (Point<2> (0.6, 0.5)); + points_glob.push_back (Point<2> (0.5, 1.0)); + points_glob.push_back (Point<2> (0.0, 1.0)); + + // Prepare cell data + std::vector > cells (1); + + switch (elm) + { + case 0: + cells[0].vertices[0] = 0; + cells[0].vertices[1] = 1; + cells[0].vertices[2] = 4; + cells[0].vertices[3] = 2; + cells[0].material_id = 0; + break; + case 1: + cells[0].vertices[0] = 4; + cells[0].vertices[1] = 2; + cells[0].vertices[2] = 5; + cells[0].vertices[3] = 3; + cells[0].material_id = 0; + break; + case 2: + cells[0].vertices[0] = 0; + cells[0].vertices[1] = 4; + cells[0].vertices[2] = 6; + cells[0].vertices[3] = 5; + cells[0].material_id = 0; + break; + } + + for (unsigned int i = 0; i < 4; ++i) + { + points.push_back (points_glob[cells[0].vertices[i]]); + cells[0].vertices[i] = i; + } + + tria.create_triangulation (points, cells, SubCellData()); + + return (0); +} + + +int plot_shapes (DoFHandler<2> &dof_handler) +{ + Vector solution (dof_handler.n_dofs ()); + std::set face_dofs; + + // Create set of all DoFs, which are on the boundary. + DoFHandler<2>::active_face_iterator face = dof_handler.begin_active_face(), + endf = dof_handler.end_face(); + std::vector face_dof_indices (dof_handler.get_fe().dofs_per_face); + for (; face!=endf; ++face) + { + face->get_dof_indices (face_dof_indices); + for (unsigned int i=0; i::iterator face_dof_iter; + for (face_dof_iter = face_dofs.begin (); + face_dof_iter != face_dofs.end (); ++face_dof_iter) + { + unsigned int dof = *face_dof_iter; + + deallog << dof << std::endl; + + solution (dof) = 1.0; + + // Test the core functionality2 + 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 (4, 0); + + data_out->write_gnuplot (deallog.get_file_stream()); + + delete data_out; + + solution (dof) = 0.0; + } + + return (0); +} + + +int main (int /*argc*/, char **/*argv*/) +{ + std::ofstream logfile (logname); + logfile.precision (PRECISION); + logfile.setf(std::ios::fixed); + deallog.attach(logfile); + deallog.depth_console(0); + deallog.threshold_double(1.e-10); + + test (); +} + + + +void check (const FiniteElement<2> &fe) +{ + Triangulation<2> tria_test; + DoFHandler<2> *dof_handler; + + deallog << "Dofs/cell " << fe.dofs_per_cell + << "Dofs/face " << fe.dofs_per_face << std::endl; + + for (unsigned int elm = 0; elm < 3; ++elm) + { + create_tria (elm, tria_test); + + // Create a DoFHandler + dof_handler = new DoFHandler<2> (tria_test); + dof_handler->distribute_dofs (fe); + + deallog << "Dofs total " << dof_handler->n_dofs () << std::endl; + + // Alloc some DoFs + Vector solution; + solution.reinit (dof_handler->n_dofs ()); + + // Project solution onto FE field + ConstraintMatrix hn_constraints; + hn_constraints.clear (); + DoFTools::make_hanging_node_constraints (*dof_handler, + hn_constraints); + hn_constraints.close (); + + MappingQ1<2> map_default; + + project (map_default, *dof_handler, hn_constraints, + QGauss6<2> (), TestMap1<2>(2), + solution); + + // 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 (4); + + data_out->write_gnuplot (deallog.get_file_stream()); + + // Now write face DoFs .. + for (unsigned int i = 0; i < 4; ++i) + { + deallog << solution(i) << std::endl; + } + + // Clean up ... + delete data_out; + delete (dof_handler); + tria_test.clear (); + } +} diff --git a/tests/fe/fe_sys_projection_01.cc b/tests/fe/fe_sys_projection_01.cc index 229dcb4451..107787871f 100644 --- a/tests/fe/fe_sys_projection_01.cc +++ b/tests/fe/fe_sys_projection_01.cc @@ -1,5 +1,5 @@ -//---------------------------- fe_sys_projection_01.cc --------------------------- -// fe_sys_projection_01.cc,v 1.3 2003/06/09 16:00:38 wolf Exp +//---------------------------- abfprojection_01.cc --------------------------- +// abfprojection_01.cc,v 1.3 2003/06/09 16:00:38 wolf Exp // Version: // // Copyright (C) 2003, 2005, 2006 by the deal.II authors @@ -9,7 +9,7 @@ // to the file deal.II/doc/license.html for the text and // further information on this license. // -//---------------------------- fe_sys_projection_01.cc --------------------------- +//---------------------------- abfprojection_01.cc --------------------------- /* * Project the function [1,1] onto a deformed grid and see whether the @@ -19,739 +19,12 @@ -#include "../tests.h" -#include +char logname[] = "abf_projection_01/output"; +#include "deformed_projection.h" -#define PRECISION 2 -#include - -std::ofstream logfile ("fe_sys_projection_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 - - -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 &, - const unsigned int ) const -{ - return (1.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)); - - // Parabolic inflow profile - for (unsigned int iCount = 0; iCount < this->n_components; iCount++) - return_value (iCount) = value (p, iCount); -} - - -/* - * Check the value of the derivative field. - */ - -void EvaluateDerivative (DoFHandler<2> *dof_handler, - Vector &solution) -{ - // This quadrature rule determines the points, where the - // derivative will be evaluated. - QGauss<2> quad (3); - FEValues<2> fe_values (dof_handler->get_fe (), quad, - UpdateFlags(update_values | - update_q_points | - update_gradients | - update_JxW_values)); - - const unsigned int n_q_points = quad.n_quadrature_points; - const unsigned int n_components = dof_handler->get_fe().n_components(); - const unsigned int dofs_per_cell = dof_handler->get_fe().dofs_per_cell; - - // Cell iterators - DoFHandler<2>::active_cell_iterator cell = dof_handler->begin_active(), - endc = dof_handler->end(); - - double err_l2 = 0, - err_hdiv = 0; - - std::vector local_dof_indices (dofs_per_cell); - - for (; cell!=endc; ++cell) - { - cell->get_dof_indices (local_dof_indices); - - fe_values.reinit (cell); - - // Get function values - std::vector > this_value(n_q_points, - Vector(n_components)); - fe_values.get_function_values (solution, this_value); - - // Get values from solution vector (For Trap.Rule) - std::vector > > - grads_here (n_q_points, std::vector > (n_components)); - fe_values.get_function_grads (solution, grads_here); - - for (unsigned int q_point=0; q_point -void create_mass_matrix (const Mapping &mapping, - const DoFHandler &dof, - const Quadrature &q, - SparseMatrix &matrix, - const Function &rhs_function, - Vector &rhs_vector, - const Function * const coefficient = 0) +void test () { - UpdateFlags update_flags = UpdateFlags(update_values | update_JxW_values | update_q_points); - if (coefficient != 0) - update_flags = UpdateFlags (update_flags | update_q_points); - - FEValues fe_values (mapping, dof.get_fe(), q, update_flags); - - const unsigned int dofs_per_cell = fe_values.dofs_per_cell, - n_q_points = fe_values.n_quadrature_points; - const FiniteElement &fe = fe_values.get_fe(); - const unsigned int n_components = fe.n_components(); - - Assert(coefficient == 0 || - coefficient->n_components==1 || - coefficient->n_components==n_components, ExcInternalError()); - - FullMatrix cell_matrix (dofs_per_cell, dofs_per_cell); - Vector cell_vector (dofs_per_cell); - std::vector coefficient_values (n_q_points); - std::vector > coefficient_vector_values (n_q_points, - Vector (n_components)); - - std::vector dof_indices (dofs_per_cell); - - std::vector > rhs_values(n_q_points, Vector(n_components)); - - typename DoFHandler::active_cell_iterator cell = dof.begin_active (), - endc = dof.end (); - for (; cell!=endc; ++cell) - { - fe_values.reinit (cell); - - cell_matrix = 0; - cell->get_dof_indices (dof_indices); - - const std::vector &weights = fe_values.get_JxW_values (); - rhs_function.vector_value_list (fe_values.get_quadrature_points(), rhs_values); - cell_vector = 0; - - if (coefficient != 0) - { - if (coefficient->n_components==1) - { - coefficient->value_list (fe_values.get_quadrature_points(), - coefficient_values); - for (unsigned int point=0; pointvector_value_list (fe_values.get_quadrature_points(), - coefficient_vector_values); - for (unsigned int point=0; point sign_change (dofs_per_cell, 1.0); - const unsigned int dofs_per_face = fe.dofs_per_face; - std::vector face_dof_indices (dofs_per_face); - - for (unsigned int f = 0; f < 2; ++f) - { - typename DoFHandler::active_face_iterator face = cell->face (f); - if (!face->at_boundary ()) - { - unsigned int nn = cell->neighbor_of_neighbor (f); - sprintf (buf, "Face %i NeigNeig %i\n", f, nn); - deallog << buf; - if (nn > 1) - { - face->get_dof_indices (face_dof_indices); - for (unsigned int j = 0; j < dofs_per_face; ++j) - { - sign_change[face_dof_indices[j]] = -1.0; - sprintf (buf, "DoF %i\n", face_dof_indices[j]); - deallog << buf; - } - } - } - } - - for (unsigned int point=0; point > val_vector (dofs_per_cell, - Vector (n_components)); - - // Precompute the component values - for (unsigned int i=0; i < dofs_per_cell; ++i) - for (unsigned int comp_i = 0; comp_i < fe.n_components (); - ++comp_i) - { - val_vector[i](comp_i) = sign_change[i] * - fe_values.shape_value_component(i,point,comp_i); - } - // Now eventually switch the sign of some of the ansatzfunctions. - // TODO - - for (unsigned int i=0; i -void create_right_hand_side (const Mapping &mapping, - const DoFHandler &dof_handler, - const Quadrature &quadrature, - const Function &rhs_function, - Vector &rhs_vector) -{ - const FiniteElement &fe = dof_handler.get_fe(); - Assert (fe.n_components() == rhs_function.n_components, - ExcInternalError()); - Assert (rhs_vector.size() == dof_handler.n_dofs(), - ExcDimensionMismatch(rhs_vector.size(), dof_handler.n_dofs())); - rhs_vector = 0; - - UpdateFlags update_flags = UpdateFlags(update_values | - update_q_points | - update_JxW_values); - FEValues fe_values (mapping, fe, quadrature, update_flags); - - const unsigned int dofs_per_cell = fe_values.dofs_per_cell, - n_q_points = fe_values.n_quadrature_points, - n_components = fe.n_components(); - - std::vector dofs (dofs_per_cell); - Vector cell_vector (dofs_per_cell); - - typename DoFHandler::active_cell_iterator cell = dof_handler.begin_active(), - endc = dof_handler.end(); - - if (n_components==1) - { - std::vector rhs_values(n_q_points); - - for (; cell!=endc; ++cell) - { - fe_values.reinit(cell); - - const std::vector &weights = fe_values.get_JxW_values (); - rhs_function.value_list (fe_values.get_quadrature_points(), rhs_values); - - cell_vector = 0; - for (unsigned int point=0; pointget_dof_indices (dofs); - - for (unsigned int i=0; i > rhs_values(n_q_points, Vector(n_components)); - - for (; cell!=endc; ++cell) - { - fe_values.reinit(cell); - - const std::vector &weights = fe_values.get_JxW_values (); - rhs_function.vector_value_list (fe_values.get_quadrature_points(), rhs_values); - - cell_vector = 0; - for (unsigned int point=0; pointget_dof_indices (dofs); - - for (unsigned int i=0; i -void project (const Mapping &mapping, - const DoFHandler &dof, - const ConstraintMatrix &constraints, - const Quadrature &quadrature, - const Function &function, - Vector &vec, - const bool enforce_zero_boundary = false, - const Quadrature & = QGauss2(), - const bool project_to_boundary_first = false) -{ - Assert (dof.get_fe().n_components() == function.n_components, - ExcInternalError()); - - const FiniteElement &fe = dof.get_fe(); - - // make up boundary values - std::map boundary_values; - - if (enforce_zero_boundary == true) - // no need to project boundary - // values, but enforce - // homogeneous boundary values - // anyway - { - // loop over all boundary faces - // to get all dof indices of - // dofs on the boundary. note - // that in 3d there are cases - // where a face is not at the - // boundary, yet one of its - // lines is, and we should - // consider the degrees of - // freedom on it as boundary - // nodes. likewise, in 2d and - // 3d there are cases where a - // cell is only at the boundary - // by one vertex. nevertheless, - // since we do not support - // boundaries with dimension - // less or equal to dim-2, each - // such boundary dof is also - // found from some other face - // that is actually wholly on - // the boundary, not only by - // one line or one vertex - typename DoFHandler::active_face_iterator face = dof.begin_active_face(), - endf = dof.end_face(); - std::vector face_dof_indices (fe.dofs_per_face); - for (; face!=endf; ++face) - if (face->at_boundary()) - { - face->get_dof_indices (face_dof_indices); - for (unsigned int i=0; i::type boundary_functions; - for (unsigned char c=0; c<255; ++c) - boundary_functions[c] = &function; - project_boundary_values (dof, boundary_functions, q_boundary, - boundary_values); -*/ - } - - - // set up mass matrix and right hand side - vec.reinit (dof.n_dofs()); - SparsityPattern sparsity(dof.n_dofs(), - dof.n_dofs(), - dof.max_couplings_between_dofs()); - DoFTools::make_sparsity_pattern (dof, sparsity); - constraints.condense (sparsity); - - SparseMatrix mass_matrix (sparsity); - Vector tmp (mass_matrix.n()); - - create_mass_matrix (mapping, dof, quadrature, mass_matrix, function, tmp); - sprintf (buf, "MM created\n"); - deallog << buf; - create_right_hand_side (mapping, dof, quadrature, function, tmp); - sprintf (buf, "RHS created\n"); - deallog << buf; - - constraints.condense (mass_matrix); - constraints.condense (tmp); - if (boundary_values.size() != 0) - MatrixTools::apply_boundary_values (boundary_values, - mass_matrix, vec, tmp, - true); - - SolverControl control(1000,1e-16); - PrimitiveVectorMemory<> memory; - SolverCG<> cg(control,memory); - - PreconditionSSOR<> prec; - prec.initialize(mass_matrix, 1.2); - // solve - cg.solve (mass_matrix, vec, tmp, prec); - - // distribute solution - constraints.distribute (vec); -} - - -int create_tria (unsigned int elm, Triangulation<2> &tria) -{ - std::vector > points_glob; - std::vector > points; - - points_glob.push_back (Point<2> (0.0, 0.0)); - points_glob.push_back (Point<2> (1.0, 0.0)); - points_glob.push_back (Point<2> (1.0, 0.5)); - points_glob.push_back (Point<2> (1.0, 1.0)); - points_glob.push_back (Point<2> (0.6, 0.5)); - points_glob.push_back (Point<2> (0.5, 1.0)); - points_glob.push_back (Point<2> (0.0, 1.0)); - - // Prepare cell data - std::vector > cells (1); - - switch (elm) - { - case 0: - cells[0].vertices[0] = 0; - cells[0].vertices[1] = 1; - cells[0].vertices[2] = 4; - cells[0].vertices[3] = 2; - cells[0].material_id = 0; - break; - case 1: - cells[0].vertices[0] = 4; - cells[0].vertices[1] = 2; - cells[0].vertices[2] = 5; - cells[0].vertices[3] = 3; - cells[0].material_id = 0; - break; - case 2: - cells[0].vertices[0] = 0; - cells[0].vertices[1] = 4; - cells[0].vertices[2] = 6; - cells[0].vertices[3] = 5; - cells[0].material_id = 0; - break; - } - - for (unsigned int i = 0; i < 4; ++i) - { - points.push_back (points_glob[cells[0].vertices[i]]); - cells[0].vertices[i] = i; - } - - tria.create_triangulation (points, cells, SubCellData()); - - return (0); -} - - -int plot_shapes (DoFHandler<2> &dof_handler) -{ - Vector solution (dof_handler.n_dofs ()); - std::set face_dofs; - - // Create set of all DoFs, which are on the boundary. - DoFHandler<2>::active_face_iterator face = dof_handler.begin_active_face(), - endf = dof_handler.end_face(); - std::vector face_dof_indices (dof_handler.get_fe().dofs_per_face); - for (; face!=endf; ++face) - { - face->get_dof_indices (face_dof_indices); - for (unsigned int i=0; i::iterator face_dof_iter; - for (face_dof_iter = face_dofs.begin (); - face_dof_iter != face_dofs.end (); ++face_dof_iter) - { - unsigned int dof = *face_dof_iter; - - sprintf (buf, "%i\n", dof); - deallog << buf; - - solution (dof) = 1.0; - - // Test the core functionality2 - 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 (4, 0); - - data_out->write_gnuplot (deallog.get_file_stream()); - - delete data_out; - - solution (dof) = 0.0; - } - - return (0); -} - - -int main (int /*argc*/, char **/*argv*/) -{ - 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; - FESystem<2> fe (FE_Q<2>(3), 2); - deallog << "Dofs/cell " << fe.dofs_per_cell - << "Dofs/face " << fe.dofs_per_face << std::endl; - - for (unsigned int elm = 0; elm < 3; ++elm) - { - create_tria (elm, tria_test); - - // Create a DoFHandler - dof_handler = new DoFHandler<2> (tria_test); - dof_handler->distribute_dofs (fe); - - deallog << "Dofs total " << dof_handler->n_dofs () << std::endl; - - // Alloc some DoFs - Vector solution; - solution.reinit (dof_handler->n_dofs ()); - - // Project solution onto FE field - ConstraintMatrix hn_constraints; - hn_constraints.clear (); - DoFTools::make_hanging_node_constraints (*dof_handler, - hn_constraints); - hn_constraints.close (); - - MappingQ1<2> map_default; - - project (map_default, *dof_handler, hn_constraints, - QGauss6<2> (), TestMap1<2>(2), - solution); - - // 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 (4); - - data_out->write_gnuplot (deallog.get_file_stream()); - - // Now write face DoFs .. - for (unsigned int i = 0; i < 4; ++i) - { - sprintf (buf, "%i %e\n", i, solution (i)); - deallog << buf; - } - - // Clean up ... - delete data_out; - delete (dof_handler); - tria_test.clear (); - } - - return (0); + check (fe); } diff --git a/tests/fe/rt_projection_01.cc b/tests/fe/rt_projection_01.cc index a28673a88d..22d0e320e7 100644 --- a/tests/fe/rt_projection_01.cc +++ b/tests/fe/rt_projection_01.cc @@ -18,741 +18,12 @@ -#include "../tests.h" -#include +char logname[] = "rt_projection_01/output"; +#include "deformed_projection.h" -#define PRECISION 2 -#include - -std::ofstream logfile ("rt_projection_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 - - -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 &, - const unsigned int ) const -{ - return (1.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)); - - // Parabolic inflow profile - for (unsigned int iCount = 0; iCount < this->n_components; iCount++) - return_value (iCount) = value (p, iCount); -} - - -/* - * Check the value of the derivative field. - */ - -void EvaluateDerivative (DoFHandler<2> *dof_handler, - Vector &solution) -{ - // This quadrature rule determines the points, where the - // derivative will be evaluated. - QGauss<2> quad (3); - FEValues<2> fe_values (dof_handler->get_fe (), quad, - UpdateFlags(update_values | - update_q_points | - update_gradients | - update_JxW_values)); - - const unsigned int n_q_points = quad.n_quadrature_points; - const unsigned int n_components = dof_handler->get_fe().n_components(); - const unsigned int dofs_per_cell = dof_handler->get_fe().dofs_per_cell; - - // Cell iterators - DoFHandler<2>::active_cell_iterator cell = dof_handler->begin_active(), - endc = dof_handler->end(); - - double err_l2 = 0, - err_hdiv = 0; - - std::vector local_dof_indices (dofs_per_cell); - - for (; cell!=endc; ++cell) - { - cell->get_dof_indices (local_dof_indices); - - fe_values.reinit (cell); - - // Get function values - std::vector > this_value(n_q_points, - Vector(n_components)); - fe_values.get_function_values (solution, this_value); - - // Get values from solution vector (For Trap.Rule) - std::vector > > - grads_here (n_q_points, std::vector > (n_components)); - fe_values.get_function_grads (solution, grads_here); - - for (unsigned int q_point=0; q_point -void create_mass_matrix (const Mapping &mapping, - const DoFHandler &dof, - const Quadrature &q, - SparseMatrix &matrix, - const Function &rhs_function, - Vector &rhs_vector, - const Function * const coefficient = 0) +void test () { - UpdateFlags update_flags = UpdateFlags(update_values | update_JxW_values | update_q_points); - if (coefficient != 0) - update_flags = UpdateFlags (update_flags | update_q_points); - - FEValues fe_values (mapping, dof.get_fe(), q, update_flags); - - const unsigned int dofs_per_cell = fe_values.dofs_per_cell, - n_q_points = fe_values.n_quadrature_points; - const FiniteElement &fe = fe_values.get_fe(); - const unsigned int n_components = fe.n_components(); - - Assert(coefficient == 0 || - coefficient->n_components==1 || - coefficient->n_components==n_components, ExcInternalError()); - - FullMatrix cell_matrix (dofs_per_cell, dofs_per_cell); - Vector cell_vector (dofs_per_cell); - std::vector coefficient_values (n_q_points); - std::vector > coefficient_vector_values (n_q_points, - Vector (n_components)); - - std::vector dof_indices (dofs_per_cell); - - std::vector > rhs_values(n_q_points, Vector(n_components)); - - typename DoFHandler::active_cell_iterator cell = dof.begin_active (), - endc = dof.end (); - for (; cell!=endc; ++cell) - { - fe_values.reinit (cell); - - cell_matrix = 0; - cell->get_dof_indices (dof_indices); - - const std::vector &weights = fe_values.get_JxW_values (); - rhs_function.vector_value_list (fe_values.get_quadrature_points(), rhs_values); - cell_vector = 0; - - if (coefficient != 0) - { - if (coefficient->n_components==1) - { - coefficient->value_list (fe_values.get_quadrature_points(), - coefficient_values); - for (unsigned int point=0; pointvector_value_list (fe_values.get_quadrature_points(), - coefficient_vector_values); - for (unsigned int point=0; point sign_change (dofs_per_cell, 1.0); - const unsigned int dofs_per_face = fe.dofs_per_face; - std::vector face_dof_indices (dofs_per_face); - - for (unsigned int f = 0; f < 2; ++f) - { - typename DoFHandler::active_face_iterator face = cell->face (f); - if (!face->at_boundary ()) - { - unsigned int nn = cell->neighbor_of_neighbor (f); - sprintf (buf, "Face %i NeigNeig %i\n", f, nn); - deallog << buf; - if (nn > 1) - { - face->get_dof_indices (face_dof_indices); - for (unsigned int j = 0; j < dofs_per_face; ++j) - { - sign_change[face_dof_indices[j]] = -1.0; - sprintf (buf, "DoF %i\n", face_dof_indices[j]); - deallog << buf; - } - } - } - } - - for (unsigned int point=0; point > val_vector (dofs_per_cell, - Vector (n_components)); - - // Precompute the component values - for (unsigned int i=0; i < dofs_per_cell; ++i) - for (unsigned int comp_i = 0; comp_i < fe.n_components (); - ++comp_i) - { - val_vector[i](comp_i) = sign_change[i] * - fe_values.shape_value_component(i,point,comp_i); - } - // Now eventually switch the sign of some of the ansatzfunctions. - // TODO - - for (unsigned int i=0; i -void create_right_hand_side (const Mapping &mapping, - const DoFHandler &dof_handler, - const Quadrature &quadrature, - const Function &rhs_function, - Vector &rhs_vector) -{ - const FiniteElement &fe = dof_handler.get_fe(); - Assert (fe.n_components() == rhs_function.n_components, - ExcInternalError()); - Assert (rhs_vector.size() == dof_handler.n_dofs(), - ExcDimensionMismatch(rhs_vector.size(), dof_handler.n_dofs())); - rhs_vector = 0; - - UpdateFlags update_flags = UpdateFlags(update_values | - update_q_points | - update_JxW_values); - FEValues fe_values (mapping, fe, quadrature, update_flags); - - const unsigned int dofs_per_cell = fe_values.dofs_per_cell, - n_q_points = fe_values.n_quadrature_points, - n_components = fe.n_components(); - - std::vector dofs (dofs_per_cell); - Vector cell_vector (dofs_per_cell); - - typename DoFHandler::active_cell_iterator cell = dof_handler.begin_active(), - endc = dof_handler.end(); - - if (n_components==1) - { - std::vector rhs_values(n_q_points); - - for (; cell!=endc; ++cell) - { - fe_values.reinit(cell); - - const std::vector &weights = fe_values.get_JxW_values (); - rhs_function.value_list (fe_values.get_quadrature_points(), rhs_values); - - cell_vector = 0; - for (unsigned int point=0; pointget_dof_indices (dofs); - - for (unsigned int i=0; i > rhs_values(n_q_points, Vector(n_components)); - - for (; cell!=endc; ++cell) - { - fe_values.reinit(cell); - - const std::vector &weights = fe_values.get_JxW_values (); - rhs_function.vector_value_list (fe_values.get_quadrature_points(), rhs_values); - - cell_vector = 0; - for (unsigned int point=0; pointget_dof_indices (dofs); - - for (unsigned int i=0; i -void project (const Mapping &mapping, - const DoFHandler &dof, - const ConstraintMatrix &constraints, - const Quadrature &quadrature, - const Function &function, - Vector &vec, - const bool enforce_zero_boundary = false, - const Quadrature & = QGauss2(), - const bool project_to_boundary_first = false) -{ - Assert (dof.get_fe().n_components() == function.n_components, - ExcInternalError()); - - const FiniteElement &fe = dof.get_fe(); - - // make up boundary values - std::map boundary_values; - - if (enforce_zero_boundary == true) - // no need to project boundary - // values, but enforce - // homogeneous boundary values - // anyway - { - // loop over all boundary faces - // to get all dof indices of - // dofs on the boundary. note - // that in 3d there are cases - // where a face is not at the - // boundary, yet one of its - // lines is, and we should - // consider the degrees of - // freedom on it as boundary - // nodes. likewise, in 2d and - // 3d there are cases where a - // cell is only at the boundary - // by one vertex. nevertheless, - // since we do not support - // boundaries with dimension - // less or equal to dim-2, each - // such boundary dof is also - // found from some other face - // that is actually wholly on - // the boundary, not only by - // one line or one vertex - typename DoFHandler::active_face_iterator face = dof.begin_active_face(), - endf = dof.end_face(); - std::vector face_dof_indices (fe.dofs_per_face); - for (; face!=endf; ++face) - if (face->at_boundary()) - { - face->get_dof_indices (face_dof_indices); - for (unsigned int i=0; i::type boundary_functions; - for (unsigned char c=0; c<255; ++c) - boundary_functions[c] = &function; - project_boundary_values (dof, boundary_functions, q_boundary, - boundary_values); -*/ - } - - - // set up mass matrix and right hand side - vec.reinit (dof.n_dofs()); - SparsityPattern sparsity(dof.n_dofs(), - dof.n_dofs(), - dof.max_couplings_between_dofs()); - DoFTools::make_sparsity_pattern (dof, sparsity); - constraints.condense (sparsity); - - SparseMatrix mass_matrix (sparsity); - Vector tmp (mass_matrix.n()); - - create_mass_matrix (mapping, dof, quadrature, mass_matrix, function, tmp); - sprintf (buf, "MM created\n"); - deallog << buf; - create_right_hand_side (mapping, dof, quadrature, function, tmp); - sprintf (buf, "RHS created\n"); - deallog << buf; - - constraints.condense (mass_matrix); - constraints.condense (tmp); - if (boundary_values.size() != 0) - MatrixTools::apply_boundary_values (boundary_values, - mass_matrix, vec, tmp, - true); - - SolverControl control(1000,1e-16); - PrimitiveVectorMemory<> memory; - SolverCG<> cg(control,memory); - - PreconditionSSOR<> prec; - prec.initialize(mass_matrix, 1.2); - // solve - cg.solve (mass_matrix, vec, tmp, prec); - - // distribute solution - constraints.distribute (vec); -} - - -int create_tria (unsigned int elm, Triangulation<2> &tria) -{ - std::vector > points_glob; - std::vector > points; - - points_glob.push_back (Point<2> (0.0, 0.0)); - points_glob.push_back (Point<2> (1.0, 0.0)); - points_glob.push_back (Point<2> (1.0, 0.5)); - points_glob.push_back (Point<2> (1.0, 1.0)); - points_glob.push_back (Point<2> (0.6, 0.5)); - points_glob.push_back (Point<2> (0.5, 1.0)); - points_glob.push_back (Point<2> (0.0, 1.0)); - - // Prepare cell data - std::vector > cells (1); - - switch (elm) - { - case 0: - cells[0].vertices[0] = 0; - cells[0].vertices[1] = 1; - cells[0].vertices[2] = 4; - cells[0].vertices[3] = 2; - cells[0].material_id = 0; - break; - case 1: - cells[0].vertices[0] = 4; - cells[0].vertices[1] = 2; - cells[0].vertices[2] = 5; - cells[0].vertices[3] = 3; - cells[0].material_id = 0; - break; - case 2: - cells[0].vertices[0] = 0; - cells[0].vertices[1] = 4; - cells[0].vertices[2] = 6; - cells[0].vertices[3] = 5; - cells[0].material_id = 0; - break; - } - - for (unsigned int i = 0; i < 4; ++i) - { - points.push_back (points_glob[cells[0].vertices[i]]); - cells[0].vertices[i] = i; - } - - tria.create_triangulation (points, cells, SubCellData()); - - return (0); -} - - -int plot_shapes (DoFHandler<2> &dof_handler) -{ - Vector solution (dof_handler.n_dofs ()); - std::set face_dofs; - - // Create set of all DoFs, which are on the boundary. - DoFHandler<2>::active_face_iterator face = dof_handler.begin_active_face(), - endf = dof_handler.end_face(); - std::vector face_dof_indices (dof_handler.get_fe().dofs_per_face); - for (; face!=endf; ++face) - { - face->get_dof_indices (face_dof_indices); - for (unsigned int i=0; i::iterator face_dof_iter; - for (face_dof_iter = face_dofs.begin (); - face_dof_iter != face_dofs.end (); ++face_dof_iter) - { - unsigned int dof = *face_dof_iter; - - sprintf (buf, "%i\n", dof); - deallog << buf; - - solution (dof) = 1.0; - - // Test the core functionality2 - 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 (4, 0); - - data_out->write_gnuplot (deallog.get_file_stream()); - - delete data_out; - - solution (dof) = 0.0; - } - - return (0); -} - - -int main (int /*argc*/, char **/*argv*/) -{ - 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; - FE_RaviartThomas<2> fe (0); - //FESystem<2> fe (FE_Q<2>(3), 2); - //FE_ABF<2> fe (0); - deallog << "Dofs/cell " << fe.dofs_per_cell - << "Dofs/face " << fe.dofs_per_face << std::endl; - - for (unsigned int elm = 0; elm < 3; ++elm) - { - create_tria (elm, tria_test); - - // Create a DoFHandler - dof_handler = new DoFHandler<2> (tria_test); - dof_handler->distribute_dofs (fe); - - deallog << "Dofs total " << dof_handler->n_dofs () << std::endl; - - // Alloc some DoFs - Vector solution; - solution.reinit (dof_handler->n_dofs ()); - - // Project solution onto FE field - ConstraintMatrix hn_constraints; - hn_constraints.clear (); - DoFTools::make_hanging_node_constraints (*dof_handler, - hn_constraints); - hn_constraints.close (); - - MappingQ1<2> map_default; - - project (map_default, *dof_handler, hn_constraints, - QGauss6<2> (), TestMap1<2>(2), - solution); - - // 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 (4); - - data_out->write_gnuplot (deallog.get_file_stream()); - - // Now write face DoFs .. - for (unsigned int i = 0; i < 4; ++i) - { - sprintf (buf, "%i %e\n", i, solution (i)); - deallog << buf; - } - - // Clean up ... - delete data_out; - delete (dof_handler); - tria_test.clear (); - } - - return (0); + check (fe); } diff --git a/tests/fe/rt_projection_01/cmp/generic b/tests/fe/rt_projection_01/cmp/generic index a6f3e7bae9..f6800890d1 100644 --- a/tests/fe/rt_projection_01/cmp/generic +++ b/tests/fe/rt_projection_01/cmp/generic @@ -2,8 +2,8 @@ DEAL::Dofs/cell 4Dofs/face 1 DEAL::Dofs total 4 DEAL::MM created -RHS created -Starting value 0.75 +DEAL::RHS created +DEAL:cg::Starting value 0.75 DEAL:cg::Convergence step 4 value 0 # This file was generated by the deal.II library. @@ -43,14 +43,14 @@ DEAL:cg::Convergence step 4 value 0 1.00 0.50 1.00 1.00 -DEAL::0 -1.000000e-01 -1 5.000000e-01 -2 1.000000e+00 -3 4.000000e-01 -Dofs total 4 +DEAL::-0.10 +DEAL::0.50 +DEAL::1.00 +DEAL::0.40 +DEAL::Dofs total 4 DEAL::MM created -RHS created -Starting value 0.45 +DEAL::RHS created +DEAL:cg::Starting value 0.45 DEAL:cg::Convergence step 4 value 0 # This file was generated by the deal.II library. @@ -90,14 +90,14 @@ DEAL:cg::Convergence step 4 value 0 1.00 1.00 1.00 1.00 -DEAL::0 6.000000e-01 -1 5.000000e-01 -2 4.000000e-01 -3 5.000000e-01 -Dofs total 4 +DEAL::0.60 +DEAL::0.50 +DEAL::0.40 +DEAL::0.50 +DEAL::Dofs total 4 DEAL::MM created -RHS created -Starting value 0.75 +DEAL::RHS created +DEAL:cg::Starting value 0.75 DEAL:cg::Convergence step 4 value 0 # This file was generated by the deal.II library. @@ -137,7 +137,7 @@ DEAL:cg::Convergence step 4 value 0 0.50 1.00 1.00 1.00 -DEAL::0 1.000000e+00 -1 6.000000e-01 -2 1.000000e-01 -3 5.000000e-01 +DEAL::1.00 +DEAL::0.60 +DEAL::0.10 +DEAL::0.50