From: bangerth Date: Mon, 31 Jul 2006 21:26:56 +0000 (+0000) Subject: New test for the abf element X-Git-Url: https://gitweb.dealii.org/cgi-bin/gitweb.cgi?a=commitdiff_plain;h=8f98a66f8b448fcc5c7d1b8bdf2515eee9c4b08c;p=dealii-svn.git New test for the abf element git-svn-id: https://svn.dealii.org/trunk@13552 0785d39b-7218-0410-832d-ea1e28bc413d --- diff --git a/tests/fe/Makefile b/tests/fe/Makefile index b05e6093e0..7912ef45e6 100644 --- a/tests/fe/Makefile +++ b/tests/fe/Makefile @@ -28,6 +28,7 @@ tests_x=fe_data_test traits fe_tools fe_tools_test mapping \ q_* \ interpolate_q1 \ rt_* \ + abf_* \ rtdiff \ rtn_* \ interpolate_rt interpolate_rtn \ diff --git a/tests/fe/abf_01.cc b/tests/fe/abf_01.cc new file mode 100644 index 0000000000..2e132d096f --- /dev/null +++ b/tests/fe/abf_01.cc @@ -0,0 +1,645 @@ +//---------------------------- abf_01.cc --------------------------- +// abf_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. +// +//---------------------------- abf_01.cc --------------------------- + + +// Show the shape functions of the Raviart-Thomas element on the unit cell +// Plots are gnuplot compatible if lines with desired prefix are selected. + +#include "../tests.h" +#include +#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 + + + + + +/* + * 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_pointget_dof_indices (face_dof_indices); + for (unsigned int j = 0; j < dofs_per_face; ++j) + { + sign_change[f * dofs_per_face + j] = -1.0; + printf ("DoF %i\n", face_dof_indices[j]); + } + } + } + } + */ + + 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); + //create_right_hand_side (mapping, dof, quadrature, function, tmp); + //printf ("RHS created\n"); + + 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_alternate_unitsquare (Triangulation<2> &tria) +{ + std::vector > points; + + points.push_back (Point<2> (0.0, 0.0)); + points.push_back (Point<2> (1.0, 0.0)); + points.push_back (Point<2> (1.0, 0.5)); + points.push_back (Point<2> (1.0, 1.0)); + points.push_back (Point<2> (0.6, 0.5)); + points.push_back (Point<2> (0.5, 1.0)); + points.push_back (Point<2> (0.0, 1.0)); + + //points.push_back (Point<2> (0.0, 0.001)); + + // Prepare cell data + std::vector > cells (3); + 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; + + cells[1].vertices[0] = 4; + cells[1].vertices[1] = 2; + cells[1].vertices[2] = 5; + cells[1].vertices[3] = 3; + cells[1].material_id = 0; + + cells[2].vertices[0] = 0; + //cells[2].vertices[0] = 7; + cells[2].vertices[1] = 4; + cells[2].vertices[2] = 6; + cells[2].vertices[3] = 5; + cells[2].material_id = 0; + + tria.create_triangulation (points, cells, SubCellData()); + + return (0); +} + + + +int main (int /*argc*/, char **/*argv*/) +{ + std::ofstream logfile ("abf_01/output"); + 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; + + create_alternate_unitsquare (tria_test); + + for (Triangulation<2>::active_cell_iterator cell = tria_test.begin_active(); + cell != tria_test.end(); ++cell) + { + deallog << "Cell " << cell << std::endl; + for (unsigned int v=0; v<4; ++v) + deallog << " " << cell->vertex(v) << std::endl; + } + + +// tria_test.refine_global (1); +// tria_test.distort_random (0.25); + + FE_ABF<2> fe (0); + deallog << "Dofs/cell " << fe.dofs_per_cell + << "Dofs/face " << fe.dofs_per_face << std::endl; + + DoFHandler<2> *dof_handler; + dof_handler = new DoFHandler<2> (tria_test); + dof_handler->distribute_dofs (fe); + + deallog << "Dofs total " << dof_handler->n_dofs () << std::endl; + + Vector solution(dof_handler->n_dofs ()); + solution = 1; + + // 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> (), ConstantFunction<2>(1., 2), + solution); + + EvaluateDerivative (dof_handler, solution); + solution.print (deallog.get_file_stream()); + + 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()); + + delete data_out; + + delete (dof_handler); +}