-//---------------------------- 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
// 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
-#include "../tests.h"
-#include <base/logstream.h>
+char logname[] = "abf_projection_01/output";
+#include "deformed_projection.h"
-#define PRECISION 2
-#include <fstream>
-
-std::ofstream logfile ("abf_projection_01/output");
-
-char buf[1000];
-#include <fstream>
-
-#include <grid/grid_generator.h>
-#include <grid/grid_in.h>
-#include <grid/grid_out.h>
-#include <grid/tria_boundary_lib.h>
-
-#include <fe/mapping_q.h>
-#include <fe/mapping_q1_eulerian.h>
-
-#include <numerics/vectors.h>
-#include <numerics/matrices.h>
-#include <numerics/data_out.h>
-#include <lac/vector.h>
-#include <lac/solver_cg.h>
-#include <lac/precondition.h>
-#include <lac/vector_memory.h>
-#include <lac/sparsity_pattern.h>
-#include <lac/sparse_matrix.h>
-
-#include <base/function.h>
-#include <base/quadrature_lib.h>
-#include <dofs/dof_constraints.h>
-#include <dofs/dof_tools.h>
-
-#include <fe/fe_system.h>
-#include <fe/fe.h>
-#include <fe/fe_q.h>
-#include <fe/fe_values.h>
-#include <fe/fe_raviart_thomas.h>
-#include <fe/fe_abf.h>
-#include <fe/fe_dgq.h>
-
-
-template <int dim>
-class TestMap1 : public Function<dim>
-{
- public:
- TestMap1 (const unsigned int n_components) :
- Function<dim> (n_components)
- {};
-
- virtual ~TestMap1 () {};
-
- virtual double value (const Point<dim> &p,
- const unsigned int component = 0) const;
-
- void vector_value (const Point<dim> &p,
- Vector<double> &return_value) const;
-};
-
-
-
-template <int dim>
-double
-TestMap1<dim>::value (const Point<dim> &,
- const unsigned int ) const
-{
- return (1.0);
-}
-
-
-template <int dim>
-void TestMap1<dim>::vector_value (const Point<dim> &p,
- Vector<double> &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<double> &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<unsigned int> 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<Vector<double> > this_value(n_q_points,
- Vector<double>(n_components));
- fe_values.get_function_values (solution, this_value);
-
- // Get values from solution vector (For Trap.Rule)
- std::vector<std::vector<Tensor<1,2> > >
- grads_here (n_q_points, std::vector<Tensor<1,2> > (n_components));
- fe_values.get_function_grads (solution, grads_here);
-
- for (unsigned int q_point=0; q_point<n_q_points; ++q_point)
- {
- // double u0 = this_value[q_point](0);
- //double v0 = this_value[q_point](1);
-
- double u0 = 0;
- double v0 = 0;
- for (unsigned int i = 0; i < dofs_per_cell; ++i)
- {
- u0 += (solution (local_dof_indices[i]) *
- fe_values.shape_value_component(i, q_point, 0));
- v0 += (solution (local_dof_indices[i]) *
- fe_values.shape_value_component(i, q_point, 1));
- }
- u0 -= 1.0;
- v0 -= 1.0;
-
- double dudx = grads_here[q_point][0][0];
- double dvdy = grads_here[q_point][1][1];
-
- /*
- sprintf (buf, "QP %i %e %e %e %e %e %e Div %e\n", q_point, u0, v0,
- dudx, dudy, dvdx, dvdy, dudx + dvdy);
- */
-
- err_l2 += (u0 * u0 + v0 * v0) * fe_values.JxW (q_point);
- err_hdiv += (dudx + dvdy) * (dudx + dvdy) * fe_values.JxW (q_point);
- }
- }
-
- sprintf (buf, "L2-Err %e Hdiv-Err %e\n", pow (err_l2, 0.5),
- pow (err_hdiv, 0.5));
- deallog << buf;
-}
-
-
-template <int dim>
-void create_mass_matrix (const Mapping<dim> &mapping,
- const DoFHandler<dim> &dof,
- const Quadrature<dim> &q,
- SparseMatrix<double> &matrix,
- const Function<dim> &rhs_function,
- Vector<double> &rhs_vector,
- const Function<dim> * 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<dim> 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<dim> &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<double> cell_matrix (dofs_per_cell, dofs_per_cell);
- Vector<double> cell_vector (dofs_per_cell);
- std::vector<double> coefficient_values (n_q_points);
- std::vector<Vector<double> > coefficient_vector_values (n_q_points,
- Vector<double> (n_components));
-
- std::vector<unsigned int> dof_indices (dofs_per_cell);
-
- std::vector<Vector<double> > rhs_values(n_q_points, Vector<double>(n_components));
-
- typename DoFHandler<dim>::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<double> &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; point<n_q_points; ++point)
- {
- const double weight = fe_values.JxW(point);
- for (unsigned int i=0; i<dofs_per_cell; ++i)
- {
- const double v = fe_values.shape_value(i,point);
- for (unsigned int j=0; j<dofs_per_cell; ++j)
- {
- const double u = fe_values.shape_value(j,point);
-
- if ((n_components==1) ||
- (fe.system_to_component_index(i).first ==
- fe.system_to_component_index(j).first))
- cell_matrix(i,j) +=
- (u * v * weight * coefficient_values[point]);
- }
- }
- }
- }
- else
- {
- coefficient->vector_value_list (fe_values.get_quadrature_points(),
- coefficient_vector_values);
- for (unsigned int point=0; point<n_q_points; ++point)
- {
- const double weight = fe_values.JxW(point);
- for (unsigned int i=0; i<dofs_per_cell; ++i)
- {
- const double v = fe_values.shape_value(i,point);
- const unsigned int component_i=
- fe.system_to_component_index(i).first;
- for (unsigned int j=0; j<dofs_per_cell; ++j)
- {
- const double u = fe_values.shape_value(j,point);
- if ((n_components==1) ||
- (fe.system_to_component_index(j).first == component_i))
- cell_matrix(i,j) +=
- (u * v * weight *
- coefficient_vector_values[point](component_i));
- }
- }
- }
- }
- }
- else
- {
- // Compute eventual sign changes depending on the neighborhood
- // between two faces.
- std::vector<double> sign_change (dofs_per_cell, 1.0);
- const unsigned int dofs_per_face = fe.dofs_per_face;
- std::vector<unsigned int> face_dof_indices (dofs_per_face);
-
- for (unsigned int f = 0; f < 2; ++f)
- {
- typename DoFHandler<dim>::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<n_q_points; ++point)
- {
- const double weight = fe_values.JxW(point);
- // const double weight = q.weight(point);
-
- std::vector<Vector<double> > val_vector (dofs_per_cell,
- Vector<double> (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<dofs_per_cell; ++i)
- for (unsigned int comp_i = 0; comp_i < fe.n_components ();
- ++comp_i)
- if (fe.get_nonzero_components(i)[comp_i] == true)
- {
- const double v = val_vector[i](comp_i);
- for (unsigned int j=0; j<dofs_per_cell; ++j)
- for (unsigned int comp_j = 0;
- comp_j < fe.n_components (); ++comp_j)
- if (fe.get_nonzero_components(j)[comp_j] == true)
- {
- const double u = val_vector[j](comp_j);
- if ((n_components==1) ||
- (comp_i == comp_j))
- cell_matrix(i,j) += (u * v * weight);
- }
- }
-
-
- for (unsigned int i=0; i<dofs_per_cell; ++i)
- for (unsigned int comp_i = 0; comp_i < fe.n_components ();
- ++comp_i)
- if (fe.get_nonzero_components(i)[comp_i] == true)
- {
- cell_vector(i) += rhs_values[point](comp_i) *
- val_vector[i](comp_i) * weights[point];
- }
- }
- }
- // transfer everything into the
- // global object. lock the
- // matrix meanwhile
- for (unsigned int i=0; i<dofs_per_cell; ++i)
- for (unsigned int j=0; j<dofs_per_cell; ++j)
- if ((n_components==1) ||
- (cell_matrix (i,j) != 0.0))
-/*
- (fe.system_to_component_index(i).first ==
- fe.system_to_component_index(j).first))
-*/
- matrix.add (dof_indices[i], dof_indices[j],
- cell_matrix(i,j));
-
- for (unsigned int i=0; i<dofs_per_cell; ++i)
- rhs_vector(dof_indices[i]) += cell_vector(i);
- }
-}
-
-
-template <int dim>
-void create_right_hand_side (const Mapping<dim> &mapping,
- const DoFHandler<dim> &dof_handler,
- const Quadrature<dim> &quadrature,
- const Function<dim> &rhs_function,
- Vector<double> &rhs_vector)
-{
- const FiniteElement<dim> &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<dim> 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<unsigned int> dofs (dofs_per_cell);
- Vector<double> cell_vector (dofs_per_cell);
-
- typename DoFHandler<dim>::active_cell_iterator cell = dof_handler.begin_active(),
- endc = dof_handler.end();
-
- if (n_components==1)
- {
- std::vector<double> rhs_values(n_q_points);
-
- for (; cell!=endc; ++cell)
- {
- fe_values.reinit(cell);
-
- const std::vector<double> &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; point<n_q_points; ++point)
- for (unsigned int i=0; i<dofs_per_cell; ++i)
- cell_vector(i) += rhs_values[point] *
- fe_values.shape_value(i,point) *
- weights[point];
-
- cell->get_dof_indices (dofs);
-
- for (unsigned int i=0; i<dofs_per_cell; ++i)
- rhs_vector(dofs[i]) += cell_vector(i);
- }
- }
- else
- {
- std::vector<Vector<double> > rhs_values(n_q_points, Vector<double>(n_components));
-
- for (; cell!=endc; ++cell)
- {
- fe_values.reinit(cell);
-
- const std::vector<double> &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; point<n_q_points; ++point)
- for (unsigned int i=0; i<dofs_per_cell; ++i)
- for (unsigned int comp_i = 0; comp_i < fe.n_components ();
- ++comp_i)
- // if (fe.get_nonzero_components(i)[comp_i] == true)
- {
- cell_vector(i) += rhs_values[point](comp_i) *
- fe_values.shape_value_component(i,point,comp_i) *
- weights[point];
- }
-
- cell->get_dof_indices (dofs);
-
- for (unsigned int i=0; i<dofs_per_cell; ++i)
- rhs_vector(dofs[i]) += cell_vector(i);
- }
- }
-}
-
-
-
-//
-// This function replaces the deal.II implementation of the projection.
-// The purpose is to have more freedom in assembling the matrix.
-//
-
-template <int dim>
-void project (const Mapping<dim> &mapping,
- const DoFHandler<dim> &dof,
- const ConstraintMatrix &constraints,
- const Quadrature<dim> &quadrature,
- const Function<dim> &function,
- Vector<double> &vec,
- const bool enforce_zero_boundary = false,
- const Quadrature<dim-1> & = QGauss2<dim-1>(),
- const bool project_to_boundary_first = false)
-{
- Assert (dof.get_fe().n_components() == function.n_components,
- ExcInternalError());
-
- const FiniteElement<dim> &fe = dof.get_fe();
-
- // make up boundary values
- std::map<unsigned int,double> 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<dim>::active_face_iterator face = dof.begin_active_face(),
- endf = dof.end_face();
- std::vector<unsigned int> 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<fe.dofs_per_face; ++i)
- // enter zero boundary values
- // for all boundary nodes
- //
- // we need not care about
- // vector valued elements here,
- // since we set all components
- boundary_values[face_dof_indices[i]] = 0.;
- };
- }
- else
- // no homogeneous boundary values
- if (project_to_boundary_first == true)
- // boundary projection required
- {
-/*
- // set up a list of boundary functions for
- // the different boundary parts. We want the
- // @p{function} to hold on all parts of the
- // boundary
- typename FunctionMap<dim>::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<double> mass_matrix (sparsity);
- Vector<double> 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<Point<2> > points_glob;
- std::vector<Point<2> > 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<CellData<2> > 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<double> solution (dof_handler.n_dofs ());
- std::set<unsigned int> 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<unsigned int> 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<dof_handler.get_fe().dofs_per_face; ++i)
- face_dofs.insert (face_dof_indices[i]);
- }
-
- // Now set a 1 at the place of the different DoFs and
- // output the solution.
- std::set<unsigned int>::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<double> 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);
}
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.
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.
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.
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
--- /dev/null
+//---------------------------- 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 <base/logstream.h>
+
+#define PRECISION 2
+
+#include <fstream>
+
+#include <grid/grid_generator.h>
+#include <grid/grid_in.h>
+#include <grid/grid_out.h>
+#include <grid/tria_boundary_lib.h>
+
+#include <fe/mapping_q.h>
+#include <fe/mapping_q1_eulerian.h>
+
+#include <numerics/vectors.h>
+#include <numerics/matrices.h>
+#include <numerics/data_out.h>
+#include <lac/vector.h>
+#include <lac/solver_cg.h>
+#include <lac/precondition.h>
+#include <lac/vector_memory.h>
+#include <lac/sparsity_pattern.h>
+#include <lac/sparse_matrix.h>
+
+#include <base/function.h>
+#include <base/quadrature_lib.h>
+#include <dofs/dof_constraints.h>
+#include <dofs/dof_tools.h>
+
+#include <fe/fe_system.h>
+#include <fe/fe.h>
+#include <fe/fe_q.h>
+#include <fe/fe_values.h>
+#include <fe/fe_raviart_thomas.h>
+#include <fe/fe_abf.h>
+#include <fe/fe_dgq.h>
+
+void test();
+
+
+template <int dim>
+class TestMap1 : public Function<dim>
+{
+ public:
+ TestMap1 (const unsigned int n_components) :
+ Function<dim> (n_components)
+ {};
+
+ virtual ~TestMap1 () {};
+
+ virtual double value (const Point<dim> &p,
+ const unsigned int component = 0) const;
+
+ void vector_value (const Point<dim> &p,
+ Vector<double> &return_value) const;
+};
+
+
+
+template <int dim>
+double
+TestMap1<dim>::value (const Point<dim> &,
+ const unsigned int ) const
+{
+ return (1.0);
+}
+
+
+template <int dim>
+void TestMap1<dim>::vector_value (const Point<dim> &p,
+ Vector<double> &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<double> &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<unsigned int> 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<Vector<double> > this_value(n_q_points,
+ Vector<double>(n_components));
+ fe_values.get_function_values (solution, this_value);
+
+ // Get values from solution vector (For Trap.Rule)
+ std::vector<std::vector<Tensor<1,2> > >
+ grads_here (n_q_points, std::vector<Tensor<1,2> > (n_components));
+ fe_values.get_function_grads (solution, grads_here);
+
+ for (unsigned int q_point=0; q_point<n_q_points; ++q_point)
+ {
+ // double u0 = this_value[q_point](0);
+ //double v0 = this_value[q_point](1);
+
+ double u0 = 0;
+ double v0 = 0;
+ for (unsigned int i = 0; i < dofs_per_cell; ++i)
+ {
+ u0 += (solution (local_dof_indices[i]) *
+ fe_values.shape_value_component(i, q_point, 0));
+ v0 += (solution (local_dof_indices[i]) *
+ fe_values.shape_value_component(i, q_point, 1));
+ }
+ u0 -= 1.0;
+ v0 -= 1.0;
+
+ double dudx = grads_here[q_point][0][0];
+ double dvdy = grads_here[q_point][1][1];
+
+ /*
+ sprintf (buf, "QP %i %e %e %e %e %e %e Div %e\n", q_point, u0, v0,
+ dudx, dudy, dvdx, dvdy, dudx + dvdy);
+ */
+
+ err_l2 += (u0 * u0 + v0 * v0) * fe_values.JxW (q_point);
+ err_hdiv += (dudx + dvdy) * (dudx + dvdy) * fe_values.JxW (q_point);
+ }
+ }
+
+ deallog << "L2-Err " << pow (err_l2, 0.5)
+ << " Hdiv-Err " << pow (err_hdiv, 0.5)
+ << std::endl;
+}
+
+
+template <int dim>
+void create_mass_matrix (const Mapping<dim> &mapping,
+ const DoFHandler<dim> &dof,
+ const Quadrature<dim> &q,
+ SparseMatrix<double> &matrix,
+ const Function<dim> &rhs_function,
+ Vector<double> &rhs_vector,
+ const Function<dim> * 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<dim> 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<dim> &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<double> cell_matrix (dofs_per_cell, dofs_per_cell);
+ Vector<double> cell_vector (dofs_per_cell);
+ std::vector<double> coefficient_values (n_q_points);
+ std::vector<Vector<double> > coefficient_vector_values (n_q_points,
+ Vector<double> (n_components));
+
+ std::vector<unsigned int> dof_indices (dofs_per_cell);
+
+ std::vector<Vector<double> > rhs_values(n_q_points, Vector<double>(n_components));
+
+ typename DoFHandler<dim>::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<double> &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; point<n_q_points; ++point)
+ {
+ const double weight = fe_values.JxW(point);
+ for (unsigned int i=0; i<dofs_per_cell; ++i)
+ {
+ const double v = fe_values.shape_value(i,point);
+ for (unsigned int j=0; j<dofs_per_cell; ++j)
+ {
+ const double u = fe_values.shape_value(j,point);
+
+ if ((n_components==1) ||
+ (fe.system_to_component_index(i).first ==
+ fe.system_to_component_index(j).first))
+ cell_matrix(i,j) +=
+ (u * v * weight * coefficient_values[point]);
+ }
+ }
+ }
+ }
+ else
+ {
+ coefficient->vector_value_list (fe_values.get_quadrature_points(),
+ coefficient_vector_values);
+ for (unsigned int point=0; point<n_q_points; ++point)
+ {
+ const double weight = fe_values.JxW(point);
+ for (unsigned int i=0; i<dofs_per_cell; ++i)
+ {
+ const double v = fe_values.shape_value(i,point);
+ const unsigned int component_i=
+ fe.system_to_component_index(i).first;
+ for (unsigned int j=0; j<dofs_per_cell; ++j)
+ {
+ const double u = fe_values.shape_value(j,point);
+ if ((n_components==1) ||
+ (fe.system_to_component_index(j).first == component_i))
+ cell_matrix(i,j) +=
+ (u * v * weight *
+ coefficient_vector_values[point](component_i));
+ }
+ }
+ }
+ }
+ }
+ else
+ {
+ // Compute eventual sign changes depending on the neighborhood
+ // between two faces.
+ std::vector<double> sign_change (dofs_per_cell, 1.0);
+ const unsigned int dofs_per_face = fe.dofs_per_face;
+ std::vector<unsigned int> face_dof_indices (dofs_per_face);
+
+ for (unsigned int f = 0; f < 2; ++f)
+ {
+ typename DoFHandler<dim>::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<n_q_points; ++point)
+ {
+ const double weight = fe_values.JxW(point);
+ // const double weight = q.weight(point);
+
+ std::vector<Vector<double> > val_vector (dofs_per_cell,
+ Vector<double> (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<dofs_per_cell; ++i)
+ for (unsigned int comp_i = 0; comp_i < fe.n_components ();
+ ++comp_i)
+ if (fe.get_nonzero_components(i)[comp_i] == true)
+ {
+ const double v = val_vector[i](comp_i);
+ for (unsigned int j=0; j<dofs_per_cell; ++j)
+ for (unsigned int comp_j = 0;
+ comp_j < fe.n_components (); ++comp_j)
+ if (fe.get_nonzero_components(j)[comp_j] == true)
+ {
+ const double u = val_vector[j](comp_j);
+ if ((n_components==1) ||
+ (comp_i == comp_j))
+ cell_matrix(i,j) += (u * v * weight);
+ }
+ }
+
+
+ for (unsigned int i=0; i<dofs_per_cell; ++i)
+ for (unsigned int comp_i = 0; comp_i < fe.n_components ();
+ ++comp_i)
+ if (fe.get_nonzero_components(i)[comp_i] == true)
+ {
+ cell_vector(i) += rhs_values[point](comp_i) *
+ val_vector[i](comp_i) * weights[point];
+ }
+ }
+ }
+ // transfer everything into the
+ // global object. lock the
+ // matrix meanwhile
+ for (unsigned int i=0; i<dofs_per_cell; ++i)
+ for (unsigned int j=0; j<dofs_per_cell; ++j)
+ if ((n_components==1) ||
+ (cell_matrix (i,j) != 0.0))
+/*
+ (fe.system_to_component_index(i).first ==
+ fe.system_to_component_index(j).first))
+*/
+ matrix.add (dof_indices[i], dof_indices[j],
+ cell_matrix(i,j));
+
+ for (unsigned int i=0; i<dofs_per_cell; ++i)
+ rhs_vector(dof_indices[i]) += cell_vector(i);
+ }
+}
+
+
+template <int dim>
+void create_right_hand_side (const Mapping<dim> &mapping,
+ const DoFHandler<dim> &dof_handler,
+ const Quadrature<dim> &quadrature,
+ const Function<dim> &rhs_function,
+ Vector<double> &rhs_vector)
+{
+ const FiniteElement<dim> &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<dim> 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<unsigned int> dofs (dofs_per_cell);
+ Vector<double> cell_vector (dofs_per_cell);
+
+ typename DoFHandler<dim>::active_cell_iterator cell = dof_handler.begin_active(),
+ endc = dof_handler.end();
+
+ if (n_components==1)
+ {
+ std::vector<double> rhs_values(n_q_points);
+
+ for (; cell!=endc; ++cell)
+ {
+ fe_values.reinit(cell);
+
+ const std::vector<double> &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; point<n_q_points; ++point)
+ for (unsigned int i=0; i<dofs_per_cell; ++i)
+ cell_vector(i) += rhs_values[point] *
+ fe_values.shape_value(i,point) *
+ weights[point];
+
+ cell->get_dof_indices (dofs);
+
+ for (unsigned int i=0; i<dofs_per_cell; ++i)
+ rhs_vector(dofs[i]) += cell_vector(i);
+ }
+ }
+ else
+ {
+ std::vector<Vector<double> > rhs_values(n_q_points, Vector<double>(n_components));
+
+ for (; cell!=endc; ++cell)
+ {
+ fe_values.reinit(cell);
+
+ const std::vector<double> &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; point<n_q_points; ++point)
+ for (unsigned int i=0; i<dofs_per_cell; ++i)
+ for (unsigned int comp_i = 0; comp_i < fe.n_components ();
+ ++comp_i)
+ // if (fe.get_nonzero_components(i)[comp_i] == true)
+ {
+ cell_vector(i) += rhs_values[point](comp_i) *
+ fe_values.shape_value_component(i,point,comp_i) *
+ weights[point];
+ }
+
+ cell->get_dof_indices (dofs);
+
+ for (unsigned int i=0; i<dofs_per_cell; ++i)
+ rhs_vector(dofs[i]) += cell_vector(i);
+ }
+ }
+}
+
+
+
+//
+// This function replaces the deal.II implementation of the projection.
+// The purpose is to have more freedom in assembling the matrix.
+//
+
+template <int dim>
+void project (const Mapping<dim> &mapping,
+ const DoFHandler<dim> &dof,
+ const ConstraintMatrix &constraints,
+ const Quadrature<dim> &quadrature,
+ const Function<dim> &function,
+ Vector<double> &vec,
+ const bool enforce_zero_boundary = false,
+ const Quadrature<dim-1> & = QGauss2<dim-1>(),
+ const bool project_to_boundary_first = false)
+{
+ Assert (dof.get_fe().n_components() == function.n_components,
+ ExcInternalError());
+
+ const FiniteElement<dim> &fe = dof.get_fe();
+
+ // make up boundary values
+ std::map<unsigned int,double> 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<dim>::active_face_iterator face = dof.begin_active_face(),
+ endf = dof.end_face();
+ std::vector<unsigned int> 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<fe.dofs_per_face; ++i)
+ // enter zero boundary values
+ // for all boundary nodes
+ //
+ // we need not care about
+ // vector valued elements here,
+ // since we set all components
+ boundary_values[face_dof_indices[i]] = 0.;
+ };
+ }
+ else
+ // no homogeneous boundary values
+ if (project_to_boundary_first == true)
+ // boundary projection required
+ {
+/*
+ // set up a list of boundary functions for
+ // the different boundary parts. We want the
+ // @p{function} to hold on all parts of the
+ // boundary
+ typename FunctionMap<dim>::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<double> mass_matrix (sparsity);
+ Vector<double> 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<Point<2> > points_glob;
+ std::vector<Point<2> > 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<CellData<2> > 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<double> solution (dof_handler.n_dofs ());
+ std::set<unsigned int> 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<unsigned int> 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<dof_handler.get_fe().dofs_per_face; ++i)
+ face_dofs.insert (face_dof_indices[i]);
+ }
+
+ // Now set a 1 at the place of the different DoFs and
+ // output the solution.
+ std::set<unsigned int>::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<double> 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 ();
+ }
+}
-//---------------------------- 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
// 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
-#include "../tests.h"
-#include <base/logstream.h>
+char logname[] = "abf_projection_01/output";
+#include "deformed_projection.h"
-#define PRECISION 2
-#include <fstream>
-
-std::ofstream logfile ("fe_sys_projection_01/output");
-
-char buf[1000];
-#include <fstream>
-
-#include <grid/grid_generator.h>
-#include <grid/grid_in.h>
-#include <grid/grid_out.h>
-#include <grid/tria_boundary_lib.h>
-
-#include <fe/mapping_q.h>
-#include <fe/mapping_q1_eulerian.h>
-
-#include <numerics/vectors.h>
-#include <numerics/matrices.h>
-#include <numerics/data_out.h>
-#include <lac/vector.h>
-#include <lac/solver_cg.h>
-#include <lac/precondition.h>
-#include <lac/vector_memory.h>
-#include <lac/sparsity_pattern.h>
-#include <lac/sparse_matrix.h>
-
-#include <base/function.h>
-#include <base/quadrature_lib.h>
-#include <dofs/dof_constraints.h>
-#include <dofs/dof_tools.h>
-
-#include <fe/fe_system.h>
-#include <fe/fe.h>
-#include <fe/fe_q.h>
-#include <fe/fe_values.h>
-#include <fe/fe_raviart_thomas.h>
-#include <fe/fe_abf.h>
-#include <fe/fe_dgq.h>
-
-
-template <int dim>
-class TestMap1 : public Function<dim>
-{
- public:
- TestMap1 (const unsigned int n_components) :
- Function<dim> (n_components)
- {};
-
- virtual ~TestMap1 () {};
-
- virtual double value (const Point<dim> &p,
- const unsigned int component = 0) const;
-
- void vector_value (const Point<dim> &p,
- Vector<double> &return_value) const;
-};
-
-
-
-template <int dim>
-double
-TestMap1<dim>::value (const Point<dim> &,
- const unsigned int ) const
-{
- return (1.0);
-}
-
-
-template <int dim>
-void TestMap1<dim>::vector_value (const Point<dim> &p,
- Vector<double> &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<double> &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<unsigned int> 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<Vector<double> > this_value(n_q_points,
- Vector<double>(n_components));
- fe_values.get_function_values (solution, this_value);
-
- // Get values from solution vector (For Trap.Rule)
- std::vector<std::vector<Tensor<1,2> > >
- grads_here (n_q_points, std::vector<Tensor<1,2> > (n_components));
- fe_values.get_function_grads (solution, grads_here);
-
- for (unsigned int q_point=0; q_point<n_q_points; ++q_point)
- {
- // double u0 = this_value[q_point](0);
- //double v0 = this_value[q_point](1);
-
- double u0 = 0;
- double v0 = 0;
- for (unsigned int i = 0; i < dofs_per_cell; ++i)
- {
- u0 += (solution (local_dof_indices[i]) *
- fe_values.shape_value_component(i, q_point, 0));
- v0 += (solution (local_dof_indices[i]) *
- fe_values.shape_value_component(i, q_point, 1));
- }
- u0 -= 1.0;
- v0 -= 1.0;
-
- double dudx = grads_here[q_point][0][0];
- double dvdy = grads_here[q_point][1][1];
-
- /*
- sprintf (buf, "QP %i %e %e %e %e %e %e Div %e\n", q_point, u0, v0,
- dudx, dudy, dvdx, dvdy, dudx + dvdy);
- */
-
- err_l2 += (u0 * u0 + v0 * v0) * fe_values.JxW (q_point);
- err_hdiv += (dudx + dvdy) * (dudx + dvdy) * fe_values.JxW (q_point);
- }
- }
-
- sprintf (buf, "L2-Err %e Hdiv-Err %e\n", pow (err_l2, 0.5),
- pow (err_hdiv, 0.5));
- deallog << buf;
-}
-
-
-template <int dim>
-void create_mass_matrix (const Mapping<dim> &mapping,
- const DoFHandler<dim> &dof,
- const Quadrature<dim> &q,
- SparseMatrix<double> &matrix,
- const Function<dim> &rhs_function,
- Vector<double> &rhs_vector,
- const Function<dim> * 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<dim> 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<dim> &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<double> cell_matrix (dofs_per_cell, dofs_per_cell);
- Vector<double> cell_vector (dofs_per_cell);
- std::vector<double> coefficient_values (n_q_points);
- std::vector<Vector<double> > coefficient_vector_values (n_q_points,
- Vector<double> (n_components));
-
- std::vector<unsigned int> dof_indices (dofs_per_cell);
-
- std::vector<Vector<double> > rhs_values(n_q_points, Vector<double>(n_components));
-
- typename DoFHandler<dim>::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<double> &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; point<n_q_points; ++point)
- {
- const double weight = fe_values.JxW(point);
- for (unsigned int i=0; i<dofs_per_cell; ++i)
- {
- const double v = fe_values.shape_value(i,point);
- for (unsigned int j=0; j<dofs_per_cell; ++j)
- {
- const double u = fe_values.shape_value(j,point);
-
- if ((n_components==1) ||
- (fe.system_to_component_index(i).first ==
- fe.system_to_component_index(j).first))
- cell_matrix(i,j) +=
- (u * v * weight * coefficient_values[point]);
- }
- }
- }
- }
- else
- {
- coefficient->vector_value_list (fe_values.get_quadrature_points(),
- coefficient_vector_values);
- for (unsigned int point=0; point<n_q_points; ++point)
- {
- const double weight = fe_values.JxW(point);
- for (unsigned int i=0; i<dofs_per_cell; ++i)
- {
- const double v = fe_values.shape_value(i,point);
- const unsigned int component_i=
- fe.system_to_component_index(i).first;
- for (unsigned int j=0; j<dofs_per_cell; ++j)
- {
- const double u = fe_values.shape_value(j,point);
- if ((n_components==1) ||
- (fe.system_to_component_index(j).first == component_i))
- cell_matrix(i,j) +=
- (u * v * weight *
- coefficient_vector_values[point](component_i));
- }
- }
- }
- }
- }
- else
- {
- // Compute eventual sign changes depending on the neighborhood
- // between two faces.
- std::vector<double> sign_change (dofs_per_cell, 1.0);
- const unsigned int dofs_per_face = fe.dofs_per_face;
- std::vector<unsigned int> face_dof_indices (dofs_per_face);
-
- for (unsigned int f = 0; f < 2; ++f)
- {
- typename DoFHandler<dim>::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<n_q_points; ++point)
- {
- const double weight = fe_values.JxW(point);
- // const double weight = q.weight(point);
-
- std::vector<Vector<double> > val_vector (dofs_per_cell,
- Vector<double> (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<dofs_per_cell; ++i)
- for (unsigned int comp_i = 0; comp_i < fe.n_components ();
- ++comp_i)
- if (fe.get_nonzero_components(i)[comp_i] == true)
- {
- const double v = val_vector[i](comp_i);
- for (unsigned int j=0; j<dofs_per_cell; ++j)
- for (unsigned int comp_j = 0;
- comp_j < fe.n_components (); ++comp_j)
- if (fe.get_nonzero_components(j)[comp_j] == true)
- {
- const double u = val_vector[j](comp_j);
- if ((n_components==1) ||
- (comp_i == comp_j))
- cell_matrix(i,j) += (u * v * weight);
- }
- }
-
-
- for (unsigned int i=0; i<dofs_per_cell; ++i)
- for (unsigned int comp_i = 0; comp_i < fe.n_components ();
- ++comp_i)
- if (fe.get_nonzero_components(i)[comp_i] == true)
- {
- cell_vector(i) += rhs_values[point](comp_i) *
- val_vector[i](comp_i) * weights[point];
- }
- }
- }
- // transfer everything into the
- // global object. lock the
- // matrix meanwhile
- for (unsigned int i=0; i<dofs_per_cell; ++i)
- for (unsigned int j=0; j<dofs_per_cell; ++j)
- if ((n_components==1) ||
- (cell_matrix (i,j) != 0.0))
-/*
- (fe.system_to_component_index(i).first ==
- fe.system_to_component_index(j).first))
-*/
- matrix.add (dof_indices[i], dof_indices[j],
- cell_matrix(i,j));
-
- for (unsigned int i=0; i<dofs_per_cell; ++i)
- rhs_vector(dof_indices[i]) += cell_vector(i);
- }
-}
-
-
-template <int dim>
-void create_right_hand_side (const Mapping<dim> &mapping,
- const DoFHandler<dim> &dof_handler,
- const Quadrature<dim> &quadrature,
- const Function<dim> &rhs_function,
- Vector<double> &rhs_vector)
-{
- const FiniteElement<dim> &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<dim> 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<unsigned int> dofs (dofs_per_cell);
- Vector<double> cell_vector (dofs_per_cell);
-
- typename DoFHandler<dim>::active_cell_iterator cell = dof_handler.begin_active(),
- endc = dof_handler.end();
-
- if (n_components==1)
- {
- std::vector<double> rhs_values(n_q_points);
-
- for (; cell!=endc; ++cell)
- {
- fe_values.reinit(cell);
-
- const std::vector<double> &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; point<n_q_points; ++point)
- for (unsigned int i=0; i<dofs_per_cell; ++i)
- cell_vector(i) += rhs_values[point] *
- fe_values.shape_value(i,point) *
- weights[point];
-
- cell->get_dof_indices (dofs);
-
- for (unsigned int i=0; i<dofs_per_cell; ++i)
- rhs_vector(dofs[i]) += cell_vector(i);
- }
- }
- else
- {
- std::vector<Vector<double> > rhs_values(n_q_points, Vector<double>(n_components));
-
- for (; cell!=endc; ++cell)
- {
- fe_values.reinit(cell);
-
- const std::vector<double> &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; point<n_q_points; ++point)
- for (unsigned int i=0; i<dofs_per_cell; ++i)
- for (unsigned int comp_i = 0; comp_i < fe.n_components ();
- ++comp_i)
- // if (fe.get_nonzero_components(i)[comp_i] == true)
- {
- cell_vector(i) += rhs_values[point](comp_i) *
- fe_values.shape_value_component(i,point,comp_i) *
- weights[point];
- }
-
- cell->get_dof_indices (dofs);
-
- for (unsigned int i=0; i<dofs_per_cell; ++i)
- rhs_vector(dofs[i]) += cell_vector(i);
- }
- }
-}
-
-
-
-//
-// This function replaces the deal.II implementation of the projection.
-// The purpose is to have more freedom in assembling the matrix.
-//
-
-template <int dim>
-void project (const Mapping<dim> &mapping,
- const DoFHandler<dim> &dof,
- const ConstraintMatrix &constraints,
- const Quadrature<dim> &quadrature,
- const Function<dim> &function,
- Vector<double> &vec,
- const bool enforce_zero_boundary = false,
- const Quadrature<dim-1> & = QGauss2<dim-1>(),
- const bool project_to_boundary_first = false)
-{
- Assert (dof.get_fe().n_components() == function.n_components,
- ExcInternalError());
-
- const FiniteElement<dim> &fe = dof.get_fe();
-
- // make up boundary values
- std::map<unsigned int,double> 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<dim>::active_face_iterator face = dof.begin_active_face(),
- endf = dof.end_face();
- std::vector<unsigned int> 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<fe.dofs_per_face; ++i)
- // enter zero boundary values
- // for all boundary nodes
- //
- // we need not care about
- // vector valued elements here,
- // since we set all components
- boundary_values[face_dof_indices[i]] = 0.;
- };
- }
- else
- // no homogeneous boundary values
- if (project_to_boundary_first == true)
- // boundary projection required
- {
-/*
- // set up a list of boundary functions for
- // the different boundary parts. We want the
- // @p{function} to hold on all parts of the
- // boundary
- typename FunctionMap<dim>::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<double> mass_matrix (sparsity);
- Vector<double> 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<Point<2> > points_glob;
- std::vector<Point<2> > 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<CellData<2> > 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<double> solution (dof_handler.n_dofs ());
- std::set<unsigned int> 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<unsigned int> 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<dof_handler.get_fe().dofs_per_face; ++i)
- face_dofs.insert (face_dof_indices[i]);
- }
-
- // Now set a 1 at the place of the different DoFs and
- // output the solution.
- std::set<unsigned int>::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<double> 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);
}
-#include "../tests.h"
-#include <base/logstream.h>
+char logname[] = "rt_projection_01/output";
+#include "deformed_projection.h"
-#define PRECISION 2
-#include <fstream>
-
-std::ofstream logfile ("rt_projection_01/output");
-
-char buf[1000];
-#include <fstream>
-
-#include <grid/grid_generator.h>
-#include <grid/grid_in.h>
-#include <grid/grid_out.h>
-#include <grid/tria_boundary_lib.h>
-
-#include <fe/mapping_q.h>
-#include <fe/mapping_q1_eulerian.h>
-
-#include <numerics/vectors.h>
-#include <numerics/matrices.h>
-#include <numerics/data_out.h>
-#include <lac/vector.h>
-#include <lac/solver_cg.h>
-#include <lac/precondition.h>
-#include <lac/vector_memory.h>
-#include <lac/sparsity_pattern.h>
-#include <lac/sparse_matrix.h>
-
-#include <base/function.h>
-#include <base/quadrature_lib.h>
-#include <dofs/dof_constraints.h>
-#include <dofs/dof_tools.h>
-
-#include <fe/fe_system.h>
-#include <fe/fe.h>
-#include <fe/fe_q.h>
-#include <fe/fe_values.h>
-#include <fe/fe_raviart_thomas.h>
-#include <fe/fe_abf.h>
-#include <fe/fe_dgq.h>
-
-
-template <int dim>
-class TestMap1 : public Function<dim>
-{
- public:
- TestMap1 (const unsigned int n_components) :
- Function<dim> (n_components)
- {};
-
- virtual ~TestMap1 () {};
-
- virtual double value (const Point<dim> &p,
- const unsigned int component = 0) const;
-
- void vector_value (const Point<dim> &p,
- Vector<double> &return_value) const;
-};
-
-
-
-template <int dim>
-double
-TestMap1<dim>::value (const Point<dim> &,
- const unsigned int ) const
-{
- return (1.0);
-}
-
-
-template <int dim>
-void TestMap1<dim>::vector_value (const Point<dim> &p,
- Vector<double> &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<double> &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<unsigned int> 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<Vector<double> > this_value(n_q_points,
- Vector<double>(n_components));
- fe_values.get_function_values (solution, this_value);
-
- // Get values from solution vector (For Trap.Rule)
- std::vector<std::vector<Tensor<1,2> > >
- grads_here (n_q_points, std::vector<Tensor<1,2> > (n_components));
- fe_values.get_function_grads (solution, grads_here);
-
- for (unsigned int q_point=0; q_point<n_q_points; ++q_point)
- {
- // double u0 = this_value[q_point](0);
- //double v0 = this_value[q_point](1);
-
- double u0 = 0;
- double v0 = 0;
- for (unsigned int i = 0; i < dofs_per_cell; ++i)
- {
- u0 += (solution (local_dof_indices[i]) *
- fe_values.shape_value_component(i, q_point, 0));
- v0 += (solution (local_dof_indices[i]) *
- fe_values.shape_value_component(i, q_point, 1));
- }
- u0 -= 1.0;
- v0 -= 1.0;
-
- double dudx = grads_here[q_point][0][0];
- double dvdy = grads_here[q_point][1][1];
-
- /*
- sprintf (buf, "QP %i %e %e %e %e %e %e Div %e\n", q_point, u0, v0,
- dudx, dudy, dvdx, dvdy, dudx + dvdy);
- */
-
- err_l2 += (u0 * u0 + v0 * v0) * fe_values.JxW (q_point);
- err_hdiv += (dudx + dvdy) * (dudx + dvdy) * fe_values.JxW (q_point);
- }
- }
-
- sprintf (buf, "L2-Err %e Hdiv-Err %e\n", pow (err_l2, 0.5),
- pow (err_hdiv, 0.5));
- deallog << buf;
-}
-
-
-template <int dim>
-void create_mass_matrix (const Mapping<dim> &mapping,
- const DoFHandler<dim> &dof,
- const Quadrature<dim> &q,
- SparseMatrix<double> &matrix,
- const Function<dim> &rhs_function,
- Vector<double> &rhs_vector,
- const Function<dim> * 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<dim> 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<dim> &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<double> cell_matrix (dofs_per_cell, dofs_per_cell);
- Vector<double> cell_vector (dofs_per_cell);
- std::vector<double> coefficient_values (n_q_points);
- std::vector<Vector<double> > coefficient_vector_values (n_q_points,
- Vector<double> (n_components));
-
- std::vector<unsigned int> dof_indices (dofs_per_cell);
-
- std::vector<Vector<double> > rhs_values(n_q_points, Vector<double>(n_components));
-
- typename DoFHandler<dim>::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<double> &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; point<n_q_points; ++point)
- {
- const double weight = fe_values.JxW(point);
- for (unsigned int i=0; i<dofs_per_cell; ++i)
- {
- const double v = fe_values.shape_value(i,point);
- for (unsigned int j=0; j<dofs_per_cell; ++j)
- {
- const double u = fe_values.shape_value(j,point);
-
- if ((n_components==1) ||
- (fe.system_to_component_index(i).first ==
- fe.system_to_component_index(j).first))
- cell_matrix(i,j) +=
- (u * v * weight * coefficient_values[point]);
- }
- }
- }
- }
- else
- {
- coefficient->vector_value_list (fe_values.get_quadrature_points(),
- coefficient_vector_values);
- for (unsigned int point=0; point<n_q_points; ++point)
- {
- const double weight = fe_values.JxW(point);
- for (unsigned int i=0; i<dofs_per_cell; ++i)
- {
- const double v = fe_values.shape_value(i,point);
- const unsigned int component_i=
- fe.system_to_component_index(i).first;
- for (unsigned int j=0; j<dofs_per_cell; ++j)
- {
- const double u = fe_values.shape_value(j,point);
- if ((n_components==1) ||
- (fe.system_to_component_index(j).first == component_i))
- cell_matrix(i,j) +=
- (u * v * weight *
- coefficient_vector_values[point](component_i));
- }
- }
- }
- }
- }
- else
- {
- // Compute eventual sign changes depending on the neighborhood
- // between two faces.
- std::vector<double> sign_change (dofs_per_cell, 1.0);
- const unsigned int dofs_per_face = fe.dofs_per_face;
- std::vector<unsigned int> face_dof_indices (dofs_per_face);
-
- for (unsigned int f = 0; f < 2; ++f)
- {
- typename DoFHandler<dim>::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<n_q_points; ++point)
- {
- const double weight = fe_values.JxW(point);
- // const double weight = q.weight(point);
-
- std::vector<Vector<double> > val_vector (dofs_per_cell,
- Vector<double> (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<dofs_per_cell; ++i)
- for (unsigned int comp_i = 0; comp_i < fe.n_components ();
- ++comp_i)
- if (fe.get_nonzero_components(i)[comp_i] == true)
- {
- const double v = val_vector[i](comp_i);
- for (unsigned int j=0; j<dofs_per_cell; ++j)
- for (unsigned int comp_j = 0;
- comp_j < fe.n_components (); ++comp_j)
- if (fe.get_nonzero_components(j)[comp_j] == true)
- {
- const double u = val_vector[j](comp_j);
- if ((n_components==1) ||
- (comp_i == comp_j))
- cell_matrix(i,j) += (u * v * weight);
- }
- }
-
-
- for (unsigned int i=0; i<dofs_per_cell; ++i)
- for (unsigned int comp_i = 0; comp_i < fe.n_components ();
- ++comp_i)
- if (fe.get_nonzero_components(i)[comp_i] == true)
- {
- cell_vector(i) += rhs_values[point](comp_i) *
- val_vector[i](comp_i) * weights[point];
- }
- }
- }
- // transfer everything into the
- // global object. lock the
- // matrix meanwhile
- for (unsigned int i=0; i<dofs_per_cell; ++i)
- for (unsigned int j=0; j<dofs_per_cell; ++j)
- if ((n_components==1) ||
- (cell_matrix (i,j) != 0.0))
-/*
- (fe.system_to_component_index(i).first ==
- fe.system_to_component_index(j).first))
-*/
- matrix.add (dof_indices[i], dof_indices[j],
- cell_matrix(i,j));
-
- for (unsigned int i=0; i<dofs_per_cell; ++i)
- rhs_vector(dof_indices[i]) += cell_vector(i);
- }
-}
-
-
-template <int dim>
-void create_right_hand_side (const Mapping<dim> &mapping,
- const DoFHandler<dim> &dof_handler,
- const Quadrature<dim> &quadrature,
- const Function<dim> &rhs_function,
- Vector<double> &rhs_vector)
-{
- const FiniteElement<dim> &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<dim> 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<unsigned int> dofs (dofs_per_cell);
- Vector<double> cell_vector (dofs_per_cell);
-
- typename DoFHandler<dim>::active_cell_iterator cell = dof_handler.begin_active(),
- endc = dof_handler.end();
-
- if (n_components==1)
- {
- std::vector<double> rhs_values(n_q_points);
-
- for (; cell!=endc; ++cell)
- {
- fe_values.reinit(cell);
-
- const std::vector<double> &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; point<n_q_points; ++point)
- for (unsigned int i=0; i<dofs_per_cell; ++i)
- cell_vector(i) += rhs_values[point] *
- fe_values.shape_value(i,point) *
- weights[point];
-
- cell->get_dof_indices (dofs);
-
- for (unsigned int i=0; i<dofs_per_cell; ++i)
- rhs_vector(dofs[i]) += cell_vector(i);
- }
- }
- else
- {
- std::vector<Vector<double> > rhs_values(n_q_points, Vector<double>(n_components));
-
- for (; cell!=endc; ++cell)
- {
- fe_values.reinit(cell);
-
- const std::vector<double> &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; point<n_q_points; ++point)
- for (unsigned int i=0; i<dofs_per_cell; ++i)
- for (unsigned int comp_i = 0; comp_i < fe.n_components ();
- ++comp_i)
- // if (fe.get_nonzero_components(i)[comp_i] == true)
- {
- cell_vector(i) += rhs_values[point](comp_i) *
- fe_values.shape_value_component(i,point,comp_i) *
- weights[point];
- }
-
- cell->get_dof_indices (dofs);
-
- for (unsigned int i=0; i<dofs_per_cell; ++i)
- rhs_vector(dofs[i]) += cell_vector(i);
- }
- }
-}
-
-
-
-//
-// This function replaces the deal.II implementation of the projection.
-// The purpose is to have more freedom in assembling the matrix.
-//
-
-template <int dim>
-void project (const Mapping<dim> &mapping,
- const DoFHandler<dim> &dof,
- const ConstraintMatrix &constraints,
- const Quadrature<dim> &quadrature,
- const Function<dim> &function,
- Vector<double> &vec,
- const bool enforce_zero_boundary = false,
- const Quadrature<dim-1> & = QGauss2<dim-1>(),
- const bool project_to_boundary_first = false)
-{
- Assert (dof.get_fe().n_components() == function.n_components,
- ExcInternalError());
-
- const FiniteElement<dim> &fe = dof.get_fe();
-
- // make up boundary values
- std::map<unsigned int,double> 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<dim>::active_face_iterator face = dof.begin_active_face(),
- endf = dof.end_face();
- std::vector<unsigned int> 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<fe.dofs_per_face; ++i)
- // enter zero boundary values
- // for all boundary nodes
- //
- // we need not care about
- // vector valued elements here,
- // since we set all components
- boundary_values[face_dof_indices[i]] = 0.;
- };
- }
- else
- // no homogeneous boundary values
- if (project_to_boundary_first == true)
- // boundary projection required
- {
-/*
- // set up a list of boundary functions for
- // the different boundary parts. We want the
- // @p{function} to hold on all parts of the
- // boundary
- typename FunctionMap<dim>::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<double> mass_matrix (sparsity);
- Vector<double> 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<Point<2> > points_glob;
- std::vector<Point<2> > 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<CellData<2> > 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<double> solution (dof_handler.n_dofs ());
- std::set<unsigned int> 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<unsigned int> 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<dof_handler.get_fe().dofs_per_face; ++i)
- face_dofs.insert (face_dof_indices[i]);
- }
-
- // Now set a 1 at the place of the different DoFs and
- // output the solution.
- std::set<unsigned int>::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<double> 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);
}
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.
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.
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.
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