From 7ce437b53654da3fcd36d69fcc582d40d64bbaae Mon Sep 17 00:00:00 2001 From: wolf Date: Tue, 6 Mar 2001 10:30:36 +0000 Subject: [PATCH] Merge with mapping-branch. git-svn-id: https://svn.dealii.org/trunk@4115 0785d39b-7218-0410-832d-ea1e28bc413d --- deal.II/deal.II/include/grid/grid_out.h | 109 +++- .../include/numerics/error_estimator.h | 4 + deal.II/deal.II/include/numerics/matrices.h | 42 +- .../source/grid/grid_out.all_dimensions.cc | 57 ++- deal.II/deal.II/source/grid/grid_out.cc | 467 +++++++++++++----- deal.II/deal.II/source/numerics/assembler.cc | 8 +- deal.II/deal.II/source/numerics/base.cc | 11 +- deal.II/deal.II/source/numerics/data_out.cc | 6 +- .../deal.II/source/numerics/data_out_faces.cc | 8 +- .../source/numerics/data_out_rotation.cc | 9 +- .../deal.II/source/numerics/data_out_stack.cc | 23 +- .../numerics/derivative_approximation.cc | 8 +- .../source/numerics/error_estimator.cc | 52 +- deal.II/deal.II/source/numerics/matrices.cc | 107 ++-- deal.II/deal.II/source/numerics/vectors.cc | 296 +++++------ deal.II/examples/step-2/step-2.cc | 39 +- deal.II/examples/step-3/step-3.cc | 22 +- deal.II/examples/step-4/step-4.cc | 14 +- deal.II/examples/step-5/step-5.cc | 5 +- deal.II/examples/step-6/step-6.cc | 35 +- 20 files changed, 820 insertions(+), 502 deletions(-) diff --git a/deal.II/deal.II/include/grid/grid_out.h b/deal.II/deal.II/include/grid/grid_out.h index b9eb33ecb9..f42be24657 100644 --- a/deal.II/deal.II/include/grid/grid_out.h +++ b/deal.II/deal.II/include/grid/grid_out.h @@ -18,6 +18,7 @@ #include template class Triangulation; +template class Mapping; /** * This class provides a means to output a triangulation to a file in different @@ -172,10 +173,21 @@ class GridOut * Default: @p{false}. */ bool write_cell_numbers; + + /** + * This is the number of + * points on a boundary face, + * that are ploted + * additionally to the + * vertices of the face. + */ + unsigned int n_boundary_face_points; + /** * Constructor. */ - GnuplotFlags (const bool write_cell_number = false); + GnuplotFlags (const bool write_cell_number = false, + const unsigned int n_boundary_face_points = 2); }; /** @@ -232,6 +244,20 @@ class GridOut * be drawn in a different color? */ bool color_lines_on_user_flag; + + /** + * This is the number of + * points on a boundary face, + * that are ploted + * additionally to the + * vertices of the face. + * + * This is used if the + * mapping used is not the + * standard @p{MappingQ1} + * mapping. + */ + unsigned int n_boundary_face_points; /** * Constructor. @@ -239,7 +265,8 @@ class GridOut EpsFlagsBase (const SizeType size_type = width, const unsigned int size = 300, const double line_width = 0.5, - bool color_lines_on_user_flag = false); + const bool color_lines_on_user_flag = false, + const unsigned int n_boundary_face_points = 2); }; @@ -252,13 +279,49 @@ class GridOut * listed in the base class */ template - struct EpsFlags : public EpsFlagsBase + struct EpsFlags; +// : public EpsFlagsBase +// { + +// /** +// * Constructor. +// */ +// EpsFlags (const SizeType size_type = width, +// const unsigned int size = 300, +// const double line_width = 0.5, +// bool color_lines_on_user_flag = false, +// const unsigned int n_boundary_face_points = 2); +// }; + +// //#if (__GNUC__==2) && (__GNUC_MINOR__ < 95) +// // template <> + + struct EpsFlags<1> : public EpsFlagsBase { + + /** + * Constructor. + */ + EpsFlags (const SizeType size_type = width, + const unsigned int size = 300, + const double line_width = 0.5, + bool color_lines_on_user_flag = false, + const unsigned int n_boundary_face_points = 2); }; -//#if (__GNUC__==2) && (__GNUC_MINOR__ < 95) -// template <> - + + struct EpsFlags<2> : public EpsFlagsBase + { + + /** + * Constructor. + */ + EpsFlags (const SizeType size_type = width, + const unsigned int size = 300, + const double line_width = 0.5, + bool color_lines_on_user_flag = false, + const unsigned int n_boundary_face_points = 2); + }; /** * Flags specific to the output of * grids in three space dimensions. @@ -292,11 +355,16 @@ class GridOut /** * Constructor. */ - EpsFlags (const double azimut_angle = 60, - const double turn_angle = 30) : - azimut_angle (azimut_angle), - turn_angle (turn_angle) {}; + EpsFlags (const SizeType size_type = width, + const unsigned int size = 300, + const double line_width = 0.5, + const bool color_lines_on_user_flag = false, + const unsigned int n_boundary_face_points = 2, + const double azimut_angle = 60, + const double turn_angle = 30); }; + + //#else //# warning Not implemented for gcc2.95 //#endif @@ -355,6 +423,12 @@ class GridOut * plots of the different levels * of grid. * + * @p{mapping} is a pointer to a + * mapping used for the + * transformation of cells at the + * boundary. If zero, then use + * standard Q1 mapping. + * * Names and values of additional * flags controlling the output * can be found in the @@ -363,7 +437,8 @@ class GridOut */ template void write_gnuplot (const Triangulation &tria, - std::ostream &out); + std::ostream &out, + const Mapping *mapping=0); /** * Write the triangulation in the @@ -444,6 +519,12 @@ class GridOut * file and can be changed there * according to need. * + * @p{mapping} is a pointer to a + * mapping used for the + * transformation of cells at the + * boundary. If zero, then use + * standard Q1 mapping. + * * Names and values of additional * flags controlling the output * can be found in the @@ -456,7 +537,8 @@ class GridOut */ template void write_eps (const Triangulation &tria, - std::ostream &out); + std::ostream &out, + const Mapping *mapping=0); /** * Write data and grid to @p{out} according @@ -467,7 +549,8 @@ class GridOut template void write (const Triangulation &tria, std::ostream &out, - const OutputFormat output_format); + const OutputFormat output_format, + const Mapping *mapping=0); /** * Set the flags to be used for output diff --git a/deal.II/deal.II/include/numerics/error_estimator.h b/deal.II/deal.II/include/numerics/error_estimator.h index 4f28bcdb4d..96f53add2f 100644 --- a/deal.II/deal.II/include/numerics/error_estimator.h +++ b/deal.II/deal.II/include/numerics/error_estimator.h @@ -1,5 +1,9 @@ //---------------------------- error_estimator.h --------------------------- +// $Id$ // Version: $Name$ +// +// Copyright (C) 1998, 1999, 2000, 2001 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 diff --git a/deal.II/deal.II/include/numerics/matrices.h b/deal.II/deal.II/include/numerics/matrices.h index c81d70f4bd..32e186486d 100644 --- a/deal.II/deal.II/include/numerics/matrices.h +++ b/deal.II/deal.II/include/numerics/matrices.h @@ -197,15 +197,7 @@ class MatrixCreator /** * Assemble the mass matrix. If no * coefficient is given, it is assumed - * to be constant one. - * - * If the coefficient is constant, it - * may be more adequate to use the - * functions assembling the mass matrix - * without quadrature. However, the - * two functions have different effects - * for finite elements composed of - * several subobjects. + * to be unity. * * See the general doc of this class * for more information. @@ -219,16 +211,8 @@ class MatrixCreator * Assemble the mass matrix and a right * hand side vector. If no * coefficient is given, it is assumed - * to be constant one. + * to be unity. * - * If the coefficient is constant, it - * may be more adequate to use the - * functions assembling the mass matrix - * without quadrature. However, the - * two functions have different effects - * for finite elements composed of - * several subobjects. - * * See the general doc of this class * for more information. */ @@ -239,28 +223,6 @@ class MatrixCreator Vector &rhs_vector, const Function *a = 0); - /** - * Create the mass matrix by exact - * evaluation without using a quadrature - * formula. - * - * No right hand side may be created using - * this function. See the general doc of - * this class for more information. - * - * It is assumed that the matrix already - * has the right size. The mass matrix - * elements are summed up to the values - * previously in the matrix, so if you want - * the pure mass matrix, you have to clear - * the matrix beforehand. - * - * See the general doc of this class - * for more information. - */ - static void create_mass_matrix (const DoFHandler &dof, - SparseMatrix &matrix); - /** * Assemble the mass matrix and a right * hand side vector along the boundary. diff --git a/deal.II/deal.II/source/grid/grid_out.all_dimensions.cc b/deal.II/deal.II/source/grid/grid_out.all_dimensions.cc index fb3a567422..812385cf80 100644 --- a/deal.II/deal.II/source/grid/grid_out.all_dimensions.cc +++ b/deal.II/deal.II/source/grid/grid_out.all_dimensions.cc @@ -22,31 +22,62 @@ GridOut::UcdFlags::UcdFlags (const bool write_preamble, {} -GridOut::GnuplotFlags::GnuplotFlags (const bool write_cell_numbers) : - write_cell_numbers (write_cell_numbers) +GridOut::GnuplotFlags::GnuplotFlags (const bool write_cell_numbers, + const unsigned int n_boundary_face_points) : + write_cell_numbers (write_cell_numbers), + n_boundary_face_points(n_boundary_face_points) {} GridOut::EpsFlagsBase::EpsFlagsBase (const SizeType size_type, const unsigned int size, const double line_width, - bool color_lines_on_user_flag) : + const bool color_lines_on_user_flag, + const unsigned int n_boundary_face_points) : size_type (size_type), size (size), line_width (line_width), - color_lines_on_user_flag(color_lines_on_user_flag) + color_lines_on_user_flag(color_lines_on_user_flag), + n_boundary_face_points(n_boundary_face_points) {} -// egcs 1.1.2 does not understand this, so I made it inlined in the -// class declaration -// -// template <> -// GridOut::EpsFlags<3>::EpsFlags (const double azimut_angle, -// const double turn_angle) : -// azimut_angle (azimut_angle), -// turn_angle (turn_angle) -// {}; + +GridOut::EpsFlags<1>::EpsFlags (const SizeType size_type, + const unsigned int size, + const double line_width, + const bool color_lines_on_user_flag, + const unsigned int n_boundary_face_points): + EpsFlagsBase(size_type, size, line_width, + color_lines_on_user_flag, + n_boundary_face_points) +{} + + +GridOut::EpsFlags<2>::EpsFlags (const SizeType size_type, + const unsigned int size, + const double line_width, + const bool color_lines_on_user_flag, + const unsigned int n_boundary_face_points): + EpsFlagsBase(size_type, size, line_width, + color_lines_on_user_flag, + n_boundary_face_points) +{} + + +GridOut::EpsFlags<3>::EpsFlags (const SizeType size_type, + const unsigned int size, + const double line_width, + const bool color_lines_on_user_flag, + const unsigned int n_boundary_face_points, + const double azimut_angle, + const double turn_angle): + EpsFlagsBase(size_type, size, line_width, + color_lines_on_user_flag, + n_boundary_face_points), + azimut_angle (azimut_angle), + turn_angle (turn_angle) +{} void GridOut::set_flags (const UcdFlags &flags) diff --git a/deal.II/deal.II/source/grid/grid_out.cc b/deal.II/deal.II/source/grid/grid_out.cc index 1eda7a75a6..0720790031 100644 --- a/deal.II/deal.II/source/grid/grid_out.cc +++ b/deal.II/deal.II/source/grid/grid_out.cc @@ -13,10 +13,12 @@ #include +#include #include #include #include #include +#include #include #include @@ -222,16 +224,68 @@ void GridOut::write_ucd_faces (const Triangulation &tria, }; }; +#if deal_II_dimension==1 + +template <> +void GridOut::write_gnuplot (const Triangulation<1> &tria, + std::ostream &out, + const Mapping<1> *) +{ + AssertThrow (out, ExcIO()); + + Triangulation<1>::active_cell_iterator cell=tria.begin_active(); + const Triangulation<1>::active_cell_iterator endc=tria.end(); + for (; cell!=endc; ++cell) + { + if (gnuplot_flags.write_cell_numbers) + out << "# cell " << cell << std::endl; +//TODO: plot level and material according to switches + + out << cell->vertex(0) + << ' ' << cell->level() + << ' ' << static_cast(cell->material_id()) << std::endl + << cell->vertex(1) + << ' ' << cell->level() + << ' ' << static_cast(cell->material_id()) << std::endl + << std::endl; + break; + } + AssertThrow (out, ExcIO()); +} + + +#else template void GridOut::write_gnuplot (const Triangulation &tria, - std::ostream &out) + std::ostream &out, + const Mapping *mapping) { AssertThrow (out, ExcIO()); + const unsigned int n_points= + gnuplot_flags.n_boundary_face_points; + typename Triangulation::active_cell_iterator cell=tria.begin_active(); const typename Triangulation::active_cell_iterator endc=tria.end(); + + // if we are to treat curved + // boundaries, then generate a + // quadrature formula which will be + // used to probe boundary points at + // curved faces + QProjector *q_projector=0; + if (mapping!=0) + { + typename std::vector > boundary_points(n_points); + for (unsigned int i=0; i quadrature(boundary_points, std::vector (n_points, 1)); + q_projector = new QProjector (quadrature, false); + } + for (; cell!=endc; ++cell) { if (gnuplot_flags.write_cell_numbers) @@ -240,16 +294,24 @@ void GridOut::write_gnuplot (const Triangulation &tria, switch (dim) { case 1: - out << cell->vertex(0) - << ' ' << cell->level() - << ' ' << static_cast(cell->material_id()) << std::endl - << cell->vertex(1) - << ' ' << cell->level() - << ' ' << static_cast(cell->material_id()) << std::endl - << std::endl; - break; - + { + Assert(false, ExcInternalError()); + break; + }; + case 2: + { + if (mapping==0 || !cell->at_boundary()) + { + // write out the four + // sides of this cell + // by putting the + // four points (+ the + // initial point + // again) in a row + // and lifting the + // drawing pencil at + // the end out << cell->vertex(0) << ' ' << cell->level() << ' ' << static_cast(cell->material_id()) << std::endl @@ -267,81 +329,153 @@ void GridOut::write_gnuplot (const Triangulation &tria, << ' ' << static_cast(cell->material_id()) << std::endl << std::endl // double new line for gnuplot 3d plots << std::endl; - break; - + } + else + // cell is at boundary + // and we are to treat + // curved + // boundaries. so loop + // over all faces and + // draw them as small + // pieces of lines + { + for (unsigned int face_no=0; + face_no::faces_per_cell; ++face_no) + { + const typename Triangulation::face_iterator + face = cell->face(face_no); + if (face->at_boundary()) + { + out << face->vertex(0) + << ' ' << cell->level() + << ' ' << static_cast(cell->material_id()) + << std::endl; + + const unsigned int offset=face_no*n_points; + for (unsigned int i=0; itransform_unit_to_real_cell + (cell, q_projector->point(offset+i))) + << ' ' << cell->level() + << ' ' << static_cast(cell->material_id()) + << std::endl; + + out << face->vertex(1) + << ' ' << cell->level() + << ' ' << static_cast(cell->material_id()) + << std::endl + << std::endl + << std::endl; + } + else + { + // if, + // however, + // the face + // is not at + // the + // boundary, + // then draw + // it as + // usual + out << face->vertex(0) + << ' ' << cell->level() + << ' ' << static_cast(cell->material_id()) + << std::endl + << face->vertex(1) + << ' ' << cell->level() + << ' ' << static_cast(cell->material_id()) + << std::endl + << std::endl + << std::endl; + }; + }; + }; + + break; + }; + case 3: - // front face - out << cell->vertex(0) - << ' ' << cell->level() - << ' ' << static_cast(cell->material_id()) << std::endl - << cell->vertex(1) - << ' ' << cell->level() - << ' ' << static_cast(cell->material_id()) << std::endl - << cell->vertex(2) - << ' ' << cell->level() - << ' ' << static_cast(cell->material_id()) << std::endl - << cell->vertex(3) - << ' ' << cell->level() - << ' ' << static_cast(cell->material_id()) << std::endl - << cell->vertex(0) - << ' ' << cell->level() - << ' ' << static_cast(cell->material_id()) << std::endl - << std::endl; - // back face - out << cell->vertex(4) - << ' ' << cell->level() - << ' ' << static_cast(cell->material_id()) << std::endl - << cell->vertex(5) - << ' ' << cell->level() - << ' ' << static_cast(cell->material_id()) << std::endl - << cell->vertex(6) - << ' ' << cell->level() - << ' ' << static_cast(cell->material_id()) << std::endl - << cell->vertex(7) - << ' ' << cell->level() - << ' ' << static_cast(cell->material_id()) << std::endl - << cell->vertex(4) - << ' ' << cell->level() - << ' ' << static_cast(cell->material_id()) << std::endl - << std::endl; - - // now for the four connecting lines - out << cell->vertex(0) - << ' ' << cell->level() - << ' ' << static_cast(cell->material_id()) << std::endl - << cell->vertex(4) - << ' ' << cell->level() - << ' ' << static_cast(cell->material_id()) << std::endl - << std::endl; - out << cell->vertex(1) - << ' ' << cell->level() - << ' ' << static_cast(cell->material_id()) << std::endl - << cell->vertex(5) - << ' ' << cell->level() - << ' ' << static_cast(cell->material_id()) << std::endl - << std::endl; - out << cell->vertex(2) - << ' ' << cell->level() - << ' ' << static_cast(cell->material_id()) << std::endl - << cell->vertex(6) - << ' ' << cell->level() - << ' ' << static_cast(cell->material_id()) << std::endl - << std::endl; - out << cell->vertex(3) - << ' ' << cell->level() - << ' ' << static_cast(cell->material_id()) << std::endl - << cell->vertex(7) - << ' ' << cell->level() - << ' ' << static_cast(cell->material_id()) << std::endl - << std::endl; - break; + { +//TODO: curved boundaries in 3d gnuplot not supported + Assert (mapping == 0, ExcNotImplemented()); + + // front face + out << cell->vertex(0) + << ' ' << cell->level() + << ' ' << static_cast(cell->material_id()) << std::endl + << cell->vertex(1) + << ' ' << cell->level() + << ' ' << static_cast(cell->material_id()) << std::endl + << cell->vertex(2) + << ' ' << cell->level() + << ' ' << static_cast(cell->material_id()) << std::endl + << cell->vertex(3) + << ' ' << cell->level() + << ' ' << static_cast(cell->material_id()) << std::endl + << cell->vertex(0) + << ' ' << cell->level() + << ' ' << static_cast(cell->material_id()) << std::endl + << std::endl; + // back face + out << cell->vertex(4) + << ' ' << cell->level() + << ' ' << static_cast(cell->material_id()) << std::endl + << cell->vertex(5) + << ' ' << cell->level() + << ' ' << static_cast(cell->material_id()) << std::endl + << cell->vertex(6) + << ' ' << cell->level() + << ' ' << static_cast(cell->material_id()) << std::endl + << cell->vertex(7) + << ' ' << cell->level() + << ' ' << static_cast(cell->material_id()) << std::endl + << cell->vertex(4) + << ' ' << cell->level() + << ' ' << static_cast(cell->material_id()) << std::endl + << std::endl; + + // now for the four connecting lines + out << cell->vertex(0) + << ' ' << cell->level() + << ' ' << static_cast(cell->material_id()) << std::endl + << cell->vertex(4) + << ' ' << cell->level() + << ' ' << static_cast(cell->material_id()) << std::endl + << std::endl; + out << cell->vertex(1) + << ' ' << cell->level() + << ' ' << static_cast(cell->material_id()) << std::endl + << cell->vertex(5) + << ' ' << cell->level() + << ' ' << static_cast(cell->material_id()) << std::endl + << std::endl; + out << cell->vertex(2) + << ' ' << cell->level() + << ' ' << static_cast(cell->material_id()) << std::endl + << cell->vertex(6) + << ' ' << cell->level() + << ' ' << static_cast(cell->material_id()) << std::endl + << std::endl; + out << cell->vertex(3) + << ' ' << cell->level() + << ' ' << static_cast(cell->material_id()) << std::endl + << cell->vertex(7) + << ' ' << cell->level() + << ' ' << static_cast(cell->material_id()) << std::endl + << std::endl; + break; + }; }; }; + + if (q_projector != 0) + delete q_projector; + AssertThrow (out, ExcIO()); }; - +#endif struct LineEntry { @@ -355,10 +489,23 @@ struct LineEntry }; +#if deal_II_dimension==1 + +template <> +void GridOut::write_eps (const Triangulation<1> &, + std::ostream &, + const Mapping<1> *) +{ + Assert(false, ExcNotImplemented()); +} + + +#else template void GridOut::write_eps (const Triangulation &tria, - std::ostream &out) + std::ostream &out, + const Mapping *mapping) { typedef std::list LineList; @@ -368,15 +515,14 @@ void GridOut::write_eps (const Triangulation &tria, // in order to avoid the recurring // distinctions between // eps_flags_1, eps_flags_2, ... - const EpsFlagsBase &eps_flags_base = (dim==1 ? - (EpsFlagsBase&)eps_flags_1 : - (dim==2 ? - (EpsFlagsBase&)eps_flags_2 : - (dim==3 ? - (EpsFlagsBase&)eps_flags_3 : - *(EpsFlagsBase*)0))); + const EpsFlagsBase &eps_flags_base = (dim==2 ? + static_cast(eps_flags_2) : + (dim==3 ? + static_cast(eps_flags_3) : + *static_cast(0))); AssertThrow (out, ExcIO()); + const unsigned int n_points = eps_flags_base.n_boundary_face_points; // make up a list of lines by which // we will construct the triangulation @@ -393,36 +539,120 @@ void GridOut::write_eps (const Triangulation &tria, switch (dim) { + case 1: + { + Assert(false, ExcInternalError()); + break; + }; + case 2: { Triangulation::active_line_iterator line =tria.begin_active_line (); Triangulation::active_line_iterator endline=tria.end_line (); - + + // first treat all interior + // lines and make up a list + // of them. if curved lines + // shall not be supported + // (i.e. no mapping is + // provided), then also treat + // all other lines for (; line!=endline; ++line) - // one would expect - // make_pair(line->vertex(0), - // line->vertex(1)) - // here, but that is not - // dimension independent, since - // vertex(i) is Point, - // but we want a Point<2>. - // in fact, whenever we're here, - // the vertex is a Point, - // but the compiler does not - // know this. hopefully, the - // compiler will optimize away - // this little kludge - line_list.push_back (LineEntry(Point<2>(line->vertex(0)(0), - line->vertex(0)(1)), - Point<2>(line->vertex(1)(0), - line->vertex(1)(1)), - line->user_flag_set())); + if (mapping==0 || !line->at_boundary()) + // one would expect + // make_pair(line->vertex(0), + // line->vertex(1)) + // here, but that is not + // dimension independent, since + // vertex(i) is Point, + // but we want a Point<2>. + // in fact, whenever we're here, + // the vertex is a Point, + // but the compiler does not + // know this. hopefully, the + // compiler will optimize away + // this little kludge + line_list.push_back (LineEntry(Point<2>(line->vertex(0)(0), + line->vertex(0)(1)), + Point<2>(line->vertex(1)(0), + line->vertex(1)(1)), + line->user_flag_set())); + + // next if we are to treat + // curved boundaries + // specially, then add lines + // to the list consisting of + // pieces of the boundary + // lines + if (mapping!=0) + { + // to do so, first + // generate a sequence of + // points on a face and + // project them onto the + // faces of a unit cell + typename std::vector > boundary_points (n_points); + + for (unsigned int i=0; i quadrature (boundary_points, + std::vector (n_points, 1)); + QProjector q_projector (quadrature, false); + + // next loop over all + // boundary faces and + // generate the info from + // them + typename Triangulation::active_cell_iterator cell=tria.begin_active (); + const typename Triangulation::active_cell_iterator end=tria.end (); + for (; cell!=end; ++cell) + for (unsigned int face_no=0; face_no::faces_per_cell; ++face_no) + { + const typename Triangulation::face_iterator + face = cell->face(face_no); + + if (face->at_boundary()) + { + Point p0_dim(face->vertex(0)); + Point<2> p0 (p0_dim(0), p0_dim(1)); + + // loop over + // all pieces + // of the line + // and generate + // line-lets + const unsigned int offset=face_no*n_points; + for (unsigned int i=0; i p1_dim (mapping->transform_unit_to_real_cell + (cell, q_projector.point(offset+i))); + const Point<2> p1 (p1_dim(0), p1_dim(1)); + + line_list.push_back (LineEntry(p0, p1, + face->user_flag_set())); + p0=p1; + } + + // generate last piece + const Point p1_dim (face->vertex(1)); + const Point<2> p1 (p1_dim(0), p1_dim(1)); + line_list.push_back (LineEntry(p0, p1, + face->user_flag_set())); + }; + }; + }; break; }; case 3: { + // curved boundary output + // presently not supported +//TODO: curved boundaries in eps for 3d + Assert (mapping == 0, ExcNotImplemented()); + Triangulation::active_line_iterator line =tria.begin_active_line (); Triangulation::active_line_iterator endline=tria.end_line (); @@ -442,14 +672,8 @@ void GridOut::write_eps (const Triangulation &tria, // is in direction of the viewer, but // I am too tired at present to fix // this -//#if !((__GNUC__==2) && (__GNUC_MINOR__==95)) const double z_angle = eps_flags_3.azimut_angle; const double turn_angle = eps_flags_3.turn_angle; -//#else -//# warning Not implemented for gcc2.95 -// const double z_angle = 60; -// const double turn_angle = 30; -//#endif const double pi = 3.1415926536; const Point view_direction(-sin(z_angle * 2.*pi / 360.) * sin(turn_angle * 2.*pi / 360.), +sin(z_angle * 2.*pi / 360.) * cos(turn_angle * 2.*pi / 360.), @@ -488,6 +712,7 @@ void GridOut::write_eps (const Triangulation &tria, }; + // find out minimum and maximum x and // y coordinates to compute offsets // and scaling factors @@ -545,9 +770,9 @@ void GridOut::write_eps (const Triangulation &tria, // lower left corner << "0 0 " // upper right corner - << static_cast( (x_max-x_min) * scale ) + << static_cast( (x_max-x_min) * scale )+1 << ' ' - << static_cast( (y_max-y_min) * scale ) + << static_cast( (y_max-y_min) * scale )+1 << std::endl; // define some abbreviations to keep @@ -582,12 +807,14 @@ void GridOut::write_eps (const Triangulation &tria, AssertThrow (out, ExcIO()); }; +#endif template void GridOut::write (const Triangulation &tria, std::ostream &out, - OutputFormat output_format) + OutputFormat output_format, + const Mapping *mapping) { switch (output_format) { @@ -596,11 +823,11 @@ void GridOut::write (const Triangulation &tria, return; case gnuplot: - write_gnuplot (tria, out); + write_gnuplot (tria, out, mapping); return; case eps: - write_eps (tria, out); + write_eps (tria, out, mapping); return; }; @@ -612,10 +839,14 @@ void GridOut::write (const Triangulation &tria, // explicit instantiations template void GridOut::write_ucd (const Triangulation &, std::ostream &); +#if deal_II_dimension != 1 template void GridOut::write_gnuplot (const Triangulation &, - std::ostream &); + std::ostream &, + const Mapping *); template void GridOut::write_eps (const Triangulation &, - std::ostream &); + std::ostream &, + const Mapping *); +#endif template void GridOut::write (const Triangulation &, std::ostream &, OutputFormat); diff --git a/deal.II/deal.II/source/numerics/assembler.cc b/deal.II/deal.II/source/numerics/assembler.cc index a1001adde3..794690cbc5 100644 --- a/deal.II/deal.II/source/numerics/assembler.cc +++ b/deal.II/deal.II/source/numerics/assembler.cc @@ -19,6 +19,9 @@ #include #include #include +#include +#include + // if necessary try to work around a bug in the IBM xlC compiler #ifdef XLC_WORK_AROUND_STD_BUG @@ -26,6 +29,9 @@ using namespace std; #endif + +static MappingQ1 mapping; + template Assembler::AssemblerData::AssemblerData (const DoFHandler &dof, const bool assemble_matrix, @@ -56,7 +62,7 @@ Assembler::Assembler (Triangulation *tria, assemble_rhs (local_data->assemble_rhs), matrix(local_data->matrix), rhs_vector(local_data->rhs_vector), - fe_values (dof_handler->get_fe(), + fe_values (mapping, dof_handler->get_fe(), local_data->quadrature, local_data->update_flags) { diff --git a/deal.II/deal.II/source/numerics/base.cc b/deal.II/deal.II/source/numerics/base.cc index 0c6e160f7d..a94cdf9920 100644 --- a/deal.II/deal.II/source/numerics/base.cc +++ b/deal.II/deal.II/source/numerics/base.cc @@ -33,11 +33,6 @@ #include #include -// if necessary try to work around a bug in the IBM xlC compiler -#ifdef XLC_WORK_AROUND_STD_BUG -using namespace std; -#endif - template ProblemBase::ProblemBase () : @@ -170,7 +165,8 @@ void ProblemBase::solve () template -void ProblemBase::fill_data (DataOut &out) const { +void ProblemBase::fill_data (DataOut &out) const +{ Assert ((tria!=0) && (dof_handler!=0), ExcNoTriaSelected()); out.clear (); @@ -181,7 +177,8 @@ void ProblemBase::fill_data (DataOut &out) const { template -std::string ProblemBase::get_solution_name () const { +std::string ProblemBase::get_solution_name () const +{ return "solution"; }; diff --git a/deal.II/deal.II/source/numerics/data_out.cc b/deal.II/deal.II/source/numerics/data_out.cc index aff5d97d6f..7902c30df4 100644 --- a/deal.II/deal.II/source/numerics/data_out.cc +++ b/deal.II/deal.II/source/numerics/data_out.cc @@ -22,6 +22,10 @@ #include #include #include +#include + +//TODO: Do this more clever +static MappingQ1 mapping; #ifdef DEAL_II_USE_MT #include @@ -275,7 +279,7 @@ void DataOut::build_some_patches (Data data) QTrapez<1> q_trapez; QIterated patch_points (q_trapez, data.n_subdivisions); - FEValues fe_patch_values(dofs->get_fe(), + FEValues fe_patch_values(mapping, dofs->get_fe(), patch_points, update_values); diff --git a/deal.II/deal.II/source/numerics/data_out_faces.cc b/deal.II/deal.II/source/numerics/data_out_faces.cc index bea1beb2a6..f29a1f763d 100644 --- a/deal.II/deal.II/source/numerics/data_out_faces.cc +++ b/deal.II/deal.II/source/numerics/data_out_faces.cc @@ -22,6 +22,11 @@ #include #include +#include + +//TODO: Do this more clever +static MappingQ1 mapping; + #ifdef DEAL_II_USE_MT #include #endif @@ -36,7 +41,8 @@ void DataOutFaces::build_some_patches (Data data) QTrapez<1> q_trapez; QIterated patch_points (q_trapez, data.n_subdivisions); - FEFaceValues fe_patch_values(dofs->get_fe(), + FEFaceValues fe_patch_values(mapping, + dofs->get_fe(), patch_points, update_values); diff --git a/deal.II/deal.II/source/numerics/data_out_rotation.cc b/deal.II/deal.II/source/numerics/data_out_rotation.cc index 2709ad609b..9142d73b4a 100644 --- a/deal.II/deal.II/source/numerics/data_out_rotation.cc +++ b/deal.II/deal.II/source/numerics/data_out_rotation.cc @@ -19,9 +19,13 @@ #include #include #include +#include #include #include +//TODO: Do this more clever +static MappingQ1 mapping; + #ifdef DEAL_II_USE_MT #include #endif @@ -35,15 +39,14 @@ using namespace std; #endif - - template void DataOutRotation::build_some_patches (Data data) { QTrapez<1> q_trapez; QIterated patch_points (q_trapez, data.n_subdivisions); - FEValues fe_patch_values(dofs->get_fe(), + FEValues fe_patch_values(mapping, + dofs->get_fe(), patch_points, update_values); diff --git a/deal.II/deal.II/source/numerics/data_out_stack.cc b/deal.II/deal.II/source/numerics/data_out_stack.cc index 0bfbd64875..6aeee12a0e 100644 --- a/deal.II/deal.II/source/numerics/data_out_stack.cc +++ b/deal.II/deal.II/source/numerics/data_out_stack.cc @@ -21,6 +21,8 @@ #include #include #include +#include + // if necessary try to work around a bug in the IBM xlC compiler @@ -29,6 +31,9 @@ using namespace std; #endif +//TODO: Comment?? Use proper mapping! +static MappingQ1 mapping; + template unsigned int @@ -77,8 +82,8 @@ void DataOutStack::attach_dof_handler (const DoFHandler &dof) template -void DataOutStack::declare_data_vector (const std::string &name, - const VectorType vector_type) +void DataOutStack::declare_data_vector (const std::string &name, + const VectorType vector_type) { std::vector names; names.push_back (name); @@ -88,7 +93,7 @@ void DataOutStack::declare_data_vector (const std::string &name, template void DataOutStack::declare_data_vector (const std::vector &names, - const VectorType vector_type) + const VectorType vector_type) { // make sure this function is // not called after some parameter @@ -129,7 +134,7 @@ void DataOutStack::declare_data_vector (const std::vector &nam template template void DataOutStack::add_data_vector (const Vector &vec, - const std::string &name) + const std::string &name) { std::vector names; names.push_back (name); @@ -165,7 +170,7 @@ void DataOutStack::add_data_vector (const Vector &vec, { data_vector->data.reinit (vec.size()); std::copy (vec.begin(), vec.end(), - data_vector->data.begin()); + data_vector->data.begin()); break; }; Assert (data_vector != dof_data.end(), @@ -218,13 +223,13 @@ void DataOutStack::build_patches (const unsigned int n_subdivisions) // cell to these points QTrapez<1> q_trapez; QIterated patch_points (q_trapez, n_subdivisions); - FEValues fe_patch_values (dof_handler->get_fe(), + FEValues fe_patch_values (mapping, dof_handler->get_fe(), patch_points, update_values); const unsigned int n_q_points = patch_points.n_quadrature_points; std::vector patch_values (n_q_points); std::vector > patch_values_system (n_q_points, - Vector(n_components)); + Vector(n_components)); // add the required number of patches DataOutBase::Patch default_patch; @@ -368,7 +373,7 @@ std::vector DataOutStack::get_dataset_names () const // explicit instantiations template class DataOutStack; template void DataOutStack::add_data_vector (const Vector &vec, - const std::string &name); + const std::string &name); template void DataOutStack::add_data_vector (const Vector &vec, - const std::string &name); + const std::string &name); diff --git a/deal.II/deal.II/source/numerics/derivative_approximation.cc b/deal.II/deal.II/source/numerics/derivative_approximation.cc index d8f400fe11..7f534b38e4 100644 --- a/deal.II/deal.II/source/numerics/derivative_approximation.cc +++ b/deal.II/deal.II/source/numerics/derivative_approximation.cc @@ -1,4 +1,4 @@ -//---------------------------- derivative_approximation.cc --------------------------- +//---------------------------- gradient_estimator.cc --------------------------- // $Id$ // Version: $Name$ // @@ -21,11 +21,15 @@ #include #include #include +#include #include #include +//TODO: Comment? Use proper mapping! +static MappingQ1 mapping; + template static T sqr (const T t) @@ -257,7 +261,7 @@ DerivativeApproximation::approximate (const DoFHandler &dof_handler, Vector &derivative_norm) { QMidpoint midpoint_rule; - FEValues fe_midpoint_value (dof_handler.get_fe(), + FEValues fe_midpoint_value (mapping, dof_handler.get_fe(), midpoint_rule, UpdateFlags(DerivativeDescription::update_flags | update_q_points)); diff --git a/deal.II/deal.II/source/numerics/error_estimator.cc b/deal.II/deal.II/source/numerics/error_estimator.cc index a026e12f07..d66742d6c1 100644 --- a/deal.II/deal.II/source/numerics/error_estimator.cc +++ b/deal.II/deal.II/source/numerics/error_estimator.cc @@ -12,25 +12,25 @@ //---------------------------- error_estimator.cc --------------------------- -#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 - // if necessary try to work around a bug in the IBM xlC compiler #ifdef XLC_WORK_AROUND_STD_BUG @@ -38,6 +38,12 @@ using namespace std; #endif + +//TODO: Comment?? Proper Mapping +static MappingQ1 mapping; + + + static inline double sqr (const double x) @@ -46,8 +52,6 @@ double sqr (const double x) }; - - #if deal_II_dimension == 1 template <> @@ -104,8 +108,7 @@ KellyErrorEstimator::Data::Data(const DoFHandler &dof_ha ExcInvalidBoundaryIndicator()); for (typename FunctionMap::const_iterator i=neumann_bc.begin(); i!=neumann_bc.end(); ++i) - Assert (i->second->n_components == n_components, ExcInvalidBoundaryFunction()); - + Assert (i->second->n_components == n_components, ExcInvalidBoundaryFunction()); // Init the size of a lot of vectors // needed in the calculations once @@ -151,7 +154,7 @@ void KellyErrorEstimator::estimate (const DoFHandler &dof_handler, const FunctionMap &neumann_bc, const Vector &solution, Vector &error, - const std::vector &component_mask, + const std::vector &component_mask, const Function *coefficients, unsigned int n_threads) { @@ -201,8 +204,8 @@ void KellyErrorEstimator<1>::estimate (const DoFHandler<1> &dof_ // if no mask given: treat all components std::vector component_mask ((component_mask_.size() == 0) ? - std::vector(n_components, true) : - component_mask_); + std::vector(n_components, true) : + component_mask_); Assert (component_mask.size() == n_components, ExcInvalidComponentMask()); Assert (count(component_mask.begin(), component_mask.end(), true) > 0, ExcInvalidComponentMask()); @@ -254,7 +257,7 @@ void KellyErrorEstimator<1>::estimate (const DoFHandler<1> &dof_ // two contributions from the two // vertices of each cell. QTrapez<1> quadrature; - FEValues<1> fe_values (dof_handler.get_fe(), quadrature, update_gradients); + FEValues<1> fe_values (mapping, dof_handler.get_fe(), quadrature, update_gradients); active_cell_iterator cell = dof_handler.begin_active(); for (unsigned int cell_index=0; cell != dof_handler.end(); ++cell, ++cell_index) { @@ -373,7 +376,8 @@ void KellyErrorEstimator::estimate_some (Data &data, // need not compute all values on the // neighbor cells, so using two objects // gives us a performance gain). - FEFaceValues fe_face_values_cell (data.dof_handler.get_fe(), + FEFaceValues fe_face_values_cell (mapping, + data.dof_handler.get_fe(), data.quadrature, UpdateFlags(update_gradients | update_JxW_values | @@ -381,10 +385,12 @@ void KellyErrorEstimator::estimate_some (Data &data, (data.coefficients != 0)) ? update_q_points : 0) | update_normal_vectors)); - FEFaceValues fe_face_values_neighbor (data.dof_handler.get_fe(), + FEFaceValues fe_face_values_neighbor (mapping, + data.dof_handler.get_fe(), data.quadrature, update_gradients); - FESubfaceValues fe_subface_values (data.dof_handler.get_fe(), + FESubfaceValues fe_subface_values (mapping, + data.dof_handler.get_fe(), data.quadrature, update_gradients); @@ -564,7 +570,7 @@ void KellyErrorEstimator::estimate (const DoFHandler &do solutions, ((component_mask.size() == 0) ? std::vector(dof_handler.get_fe().n_components(), - true) : + true) : component_mask), coefficients, n_threads, @@ -989,8 +995,8 @@ integrate_over_irregular_face (Data &data, data.face_integrals.end(), ExcInternalError()); Assert (data.face_integrals[face->child(subface_no)][0] >= 0, - ExcInternalError()); - + ExcInternalError()); + for (unsigned int n=0; nchild(subface_no)][n]; }; diff --git a/deal.II/deal.II/source/numerics/matrices.cc b/deal.II/deal.II/source/numerics/matrices.cc index ec2156ee8f..b6a0cee4b9 100644 --- a/deal.II/deal.II/source/numerics/matrices.cc +++ b/deal.II/deal.II/source/numerics/matrices.cc @@ -26,17 +26,28 @@ #include #include #include +#include #include #include #include + // if necessary try to work around a bug in the IBM xlC compiler #ifdef XLC_WORK_AROUND_STD_BUG using namespace std; #endif + + +//TODO: Comment?? Use proper mapping! +static MappingQ1 mapping; + + + +//TODO: re-create the create_mass_matrix function with 2 args + template void MatrixCreator::create_mass_matrix (const DoFHandler &dof, const Quadrature &q, @@ -96,42 +107,16 @@ void MatrixCreator::create_mass_matrix (const DoFHandler &dof, -template -void MatrixCreator::create_mass_matrix (const DoFHandler &dof, - SparseMatrix &matrix) -{ - const FiniteElement &fe = dof.get_fe(); - - const unsigned int dofs_per_cell = fe.dofs_per_cell; - - FullMatrix local_mass_matrix (dofs_per_cell, dofs_per_cell); - std::vector dofs_on_this_cell (dofs_per_cell); - - DoFHandler::active_cell_iterator cell = dof.begin_active(), - endc = dof.end(); - for (; cell!=endc; ++cell) - { - cell->get_dof_indices (dofs_on_this_cell); - fe.get_local_mass_matrix (cell, local_mass_matrix); - - for (unsigned int i=0; i -void MatrixCreator<1>::create_boundary_mass_matrix (const DoFHandler<1> &, - const Quadrature<0> &, - SparseMatrix &, - const FunctionMap &, - Vector &, - std::vector &, - const Function<1> *) +void MatrixCreator<1>::create_boundary_mass_matrix (const DoFHandler<1> &, + const Quadrature<0> &, + SparseMatrix &, + const FunctionMap &, + Vector &, + std::vector &, + const Function<1> *) { Assert (false, ExcNotImplemented()); }; @@ -140,13 +125,13 @@ void MatrixCreator<1>::create_boundary_mass_matrix (const DoFHandler<1> &, template -void MatrixCreator::create_boundary_mass_matrix (const DoFHandler &dof, - const Quadrature &q, - SparseMatrix &matrix, - const FunctionMap &rhs, - Vector &rhs_vector, - std::vector &dof_to_boundary_mapping, - const Function *a) +void MatrixCreator::create_boundary_mass_matrix (const DoFHandler &dof, + const Quadrature &q, + SparseMatrix &matrix, + const FunctionMap &rhs, + Vector &rhs_vector, + std::vector &dof_to_boundary_mapping, + const Function *a) { const FiniteElement &fe = dof.get_fe(); const unsigned int n_components = fe.n_components(); @@ -184,7 +169,7 @@ void MatrixCreator::create_boundary_mass_matrix (const DoFHandler & UpdateFlags update_flags = UpdateFlags (update_values | update_JxW_values | update_q_points); - FEFaceValues fe_values (fe, q, update_flags); + FEFaceValues fe_values (mapping, fe, q, update_flags); // two variables for the coefficient, // one for the two cases indicated in @@ -452,6 +437,8 @@ void MatrixCreator::create_laplace_matrix (const DoFHandler &dof, + +//TODO: recreate this function /* template @@ -488,6 +475,8 @@ void MatrixCreator::create_level_laplace_matrix (unsigned int level, + + template void MatrixCreator::create_laplace_matrix (const DoFHandler &dof, const Quadrature &q, @@ -525,10 +514,10 @@ template template void MatrixTools::apply_boundary_values (const std::map &boundary_values, - SparseMatrix &matrix, - Vector &solution, - Vector &right_hand_side, - const bool preserve_symmetry) + SparseMatrix &matrix, + Vector &solution, + Vector &right_hand_side, + const bool preserve_symmetry) { Assert (matrix.n() == matrix.m(), ExcDimensionsDontMatch(matrix.n(), matrix.m())); @@ -542,8 +531,8 @@ MatrixTools::apply_boundary_values (const std::map &bo return; - std::map::const_iterator dof = boundary_values.begin(), - endd = boundary_values.end(); + std::map::const_iterator dof = boundary_values.begin(), + endd = boundary_values.end(); const unsigned int n_dofs = matrix.m(); const SparsityPattern &sparsity = matrix.get_sparsity_pattern(); const unsigned int *sparsity_rowstart = sparsity.get_rowstart_indices(); @@ -695,7 +684,7 @@ MatrixTools::apply_boundary_values (const std::map &bo BlockSparseMatrix &matrix, BlockVector &solution, BlockVector &right_hand_side, - const bool preserve_symmetry) + const bool preserve_symmetry) { const unsigned int blocks = matrix.n_block_rows(); @@ -724,8 +713,8 @@ MatrixTools::apply_boundary_values (const std::map &bo return; - std::map::const_iterator dof = boundary_values.begin(), - endd = boundary_values.end(); + std::map::const_iterator dof = boundary_values.begin(), + endd = boundary_values.end(); const unsigned int n_dofs = matrix.m(); const BlockSparsityPattern & sparsity_pattern = matrix.get_sparsity_pattern(); @@ -949,17 +938,17 @@ MatrixTools::apply_boundary_values (const std::map &bo [this_sparsity.get_rowstart_indices()[row]]; else p = std::lower_bound(&this_sparsity.get_column_numbers() - [this_sparsity.get_rowstart_indices()[row]+1], - &this_sparsity.get_column_numbers() - [this_sparsity.get_rowstart_indices()[row+1]], - block_index.second); + [this_sparsity.get_rowstart_indices()[row]+1], + &this_sparsity.get_column_numbers() + [this_sparsity.get_rowstart_indices()[row+1]], + block_index.second); } else p = std::lower_bound(&this_sparsity.get_column_numbers() - [this_sparsity.get_rowstart_indices()[row]], - &this_sparsity.get_column_numbers() - [this_sparsity.get_rowstart_indices()[row+1]], - block_index.second); + [this_sparsity.get_rowstart_indices()[row]], + &this_sparsity.get_column_numbers() + [this_sparsity.get_rowstart_indices()[row+1]], + block_index.second); // check whether this line has // an entry in the regarding column @@ -1205,6 +1194,8 @@ void MassMatrix::assemble (Vector &rhs, + + template LaplaceMatrix::LaplaceMatrix (const Function * const rhs, const Function * const a) : diff --git a/deal.II/deal.II/source/numerics/vectors.cc b/deal.II/deal.II/source/numerics/vectors.cc index b99173ebc9..69f15b5e67 100644 --- a/deal.II/deal.II/source/numerics/vectors.cc +++ b/deal.II/deal.II/source/numerics/vectors.cc @@ -29,6 +29,7 @@ #include #include #include +#include #include #include @@ -42,6 +43,12 @@ using namespace std; #endif + +//TODO: Comment?? Use proper mapping! +static MappingQ1 mapping; + + + static inline double sqr (const double x) { return x*x; @@ -57,6 +64,7 @@ inline double sqr_point (const Tensor<1,dim> &p) + template void VectorTools::interpolate (const DoFHandler &dof, const Function &function, @@ -83,89 +91,10 @@ void VectorTools::interpolate (const DoFHandler &dof, // avoided to evaluate // the vectorfunction multiply at // the same point on a cell. - std::vector > unit_support_points (fe.dofs_per_cell); + std::vector > unit_support_points; fe.get_unit_support_points(unit_support_points); - - // The following works well - // if @p{dofs_per_x<=1 (x=vertex,line,cell)} - // as then - // the multiple support_points - // are placed one after another. - - // find the support points - // on a cell that - // are multiply mentioned in - // @p{unit_support_points}. - // Mark the first representative - // of each multiply mentioned - // support point by setting - // @p{true} in the boolean vector - // @p{is_representative_point}. -// vector is_representative_point(fe.dofs_per_cell, false); -// is_representative_point[0]=true; -// unsigned int n_rep_points=1; -// for (unsigned int last_rep_point=0, i=1; i dofs_on_cell (fe.dofs_per_cell); -// vector > support_points (fe.dofs_per_cell); - -// vector > rep_points (n_rep_points); -// vector > function_values_at_rep_points ( -// n_rep_points, Vector(fe.n_components())); - -// for (; cell!=endc; ++cell) -// { -// // for each cell: -// // get location of finite element -// // off-points (support_points) -// fe.get_support_points (cell, support_points); - -// // pick out the representative -// // support points -// unsigned int j=0; -// for (unsigned int i=0; iget_dof_indices (dofs_on_cell); - -// // distribute function values to the -// // whole vector -// int last_rep_point = -1; -// // it holds `is_representative_point[0]=true' -// // therefore the first @p{last_rep_point} is 0 -// // and we need to start with -// // `last_rep_point = -1' -// for (unsigned int i=0; i1}, - // i.e. it is usable also for systems - // including - // FEQ3, FEQ4, FEDG_Qx. + Assert (unit_support_points.size() != 0, + ExcNonInterpolatingFE()); // Find the support points // on a cell that @@ -222,7 +151,7 @@ void VectorTools::interpolate (const DoFHandler &dof, Assert(dof_to_rep_index_table.size()==fe.dofs_per_cell, ExcInternalError()); std::vector dofs_on_cell (fe.dofs_per_cell); - std::vector > support_points (fe.dofs_per_cell); + std::vector dummy_weights (fe.dofs_per_cell); std::vector > rep_points (n_rep_points); @@ -236,12 +165,23 @@ void VectorTools::interpolate (const DoFHandler &dof, std::vector > function_values_system (n_rep_points, Vector(fe.n_components())); + // Make a quadrature rule from support points + // to feed it into FEValues + Quadrature support_quadrature(unit_support_points, dummy_weights); + +//TODO: Higher order mapping? + // Transformed support points are computed by + // FEValues + FEValues fe_values (mapping, fe, support_quadrature, update_q_points); + for (; cell!=endc; ++cell) { // for each cell: // get location of finite element - // off-points (support_points) - fe.get_support_points (cell, support_points); + // support_points + fe_values.reinit(cell); + const std::vector >& support_points = + fe_values.get_quadrature_points(); // pick out the representative // support points @@ -289,6 +229,7 @@ void VectorTools::interpolate (const DoFHandler &dof, + template void VectorTools::interpolate (const DoFHandler &dof_1, const DoFHandler &dof_2, @@ -452,17 +393,7 @@ void VectorTools::project (const DoFHandler &dof, SparseMatrix mass_matrix (sparsity); Vector tmp (mass_matrix.n()); - // try to assemble the mass matrix by exact - // integration. if this is not supported, - // then use quadrature - try - { - MatrixCreator::create_mass_matrix (dof, mass_matrix); - } - catch (FiniteElement::ExcComputationNotUseful) - { - MatrixCreator::create_mass_matrix (dof, quadrature, mass_matrix); - }; + MatrixCreator::create_mass_matrix (dof, quadrature, mass_matrix); VectorTools::create_right_hand_side (dof, quadrature, function, tmp); @@ -504,7 +435,7 @@ void VectorTools::create_right_hand_side (const DoFHandler &dof_handler, UpdateFlags update_flags = UpdateFlags(update_values | update_q_points | update_JxW_values); - FEValues fe_values (fe, quadrature, update_flags); + 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, @@ -513,8 +444,8 @@ void VectorTools::create_right_hand_side (const DoFHandler &dof_handler, std::vector dofs (dofs_per_cell); Vector cell_vector (dofs_per_cell); - DoFHandler::active_cell_iterator cell = dof_handler.begin_active(), - endc = dof_handler.end(); + typename DoFHandler::active_cell_iterator cell = dof_handler.begin_active(), + endc = dof_handler.end(); if (n_components==1) { @@ -524,8 +455,8 @@ void VectorTools::create_right_hand_side (const DoFHandler &dof_handler, { fe_values.reinit(cell); - const FullMatrix &values = fe_values.get_shape_values (); - const std::vector &weights = fe_values.get_JxW_values (); + const FullMatrix &values = fe_values.get_shape_values (); + const std::vector &weights = fe_values.get_JxW_values (); rhs_function.value_list (fe_values.get_quadrature_points(), rhs_values); cell_vector.clear(); @@ -550,8 +481,8 @@ void VectorTools::create_right_hand_side (const DoFHandler &dof_handler, { fe_values.reinit(cell); - const FullMatrix &values = fe_values.get_shape_values (); - const std::vector &weights = fe_values.get_JxW_values (); + const FullMatrix &values = fe_values.get_shape_values (); + const std::vector &weights = fe_values.get_JxW_values (); rhs_function.vector_value_list (fe_values.get_quadrature_points(), rhs_values); cell_vector.clear(); @@ -597,8 +528,8 @@ VectorTools::interpolate_boundary_values (const DoFHandler<1> &dof, // the original value or a vector // of @p{true}s const std::vector component_mask ((component_mask_.size() == 0) ? - std::vector (fe.n_components(), true) : - component_mask_); + std::vector (fe.n_components(), true) : + component_mask_); Assert (count(component_mask.begin(), component_mask.end(), true) > 0, ExcComponentMismatch()); @@ -675,14 +606,15 @@ VectorTools::interpolate_boundary_values (const DoFHandler &dof, // the original value or a vector // of @p{true}s const std::vector component_mask ((component_mask_.size() == 0) ? - std::vector (fe.n_components(), true) : - component_mask_); + std::vector (fe.n_components(), true) : + component_mask_); Assert (count(component_mask.begin(), component_mask.end(), true) > 0, ExcComponentMismatch()); - // field to store the indices of dofs + // field to store the indices std::vector face_dofs (fe.dofs_per_face, -1); std::vector > dof_locations (face_dofs.size(), Point()); + // array to store the values of // the boundary function at the // boundary points. have to arrays @@ -691,55 +623,78 @@ VectorTools::interpolate_boundary_values (const DoFHandler &dof, // respectively std::vector dof_values_scalar (fe.dofs_per_face); std::vector > dof_values_system (fe.dofs_per_face, - Vector(fe.n_components())); - - DoFHandler::active_face_iterator face = dof.begin_active_face(), - endf = dof.end_face(); - for (; face!=endf; ++face) - if (boundary_component == face->boundary_indicator()) - // face is of the right component - { - // get indices, physical location and - // boundary values of dofs on this - // face - face->get_dof_indices (face_dofs); - fe.get_face_support_points (face, dof_locations); + Vector(fe.n_components())); - if (fe_is_system) + // next generate a quadrature rule + // on the face from the unit + // support points. this wil be used + // to obtain the quadrature points + // on the real cell's face + typename std::vector > unit_support_points; + fe.get_unit_face_support_points(unit_support_points); + // check whether there are support + // points on the face, if not, then + // this FE does not allow to be + // used in this function + Assert (unit_support_points.size() != 0, ExcNonInterpolatingFE()); + + std::vector dummy_weights (unit_support_points.size()); + Quadrature aux_quad (unit_support_points, dummy_weights); + FEFaceValues fe_values (mapping, fe, aux_quad, update_q_points); + + DoFHandler::active_cell_iterator cell = dof.begin_active(), + endc = dof.end(); + for (; cell!=endc; ++cell) + for (unsigned int face_no = 0; face_no < GeometryInfo::faces_per_cell; + ++face_no) + { + DoFHandler::face_iterator face = cell->face(face_no); + if (boundary_component == face->boundary_indicator()) + // face is of the right component { - boundary_function.vector_value_list (dof_locations, dof_values_system); - - // enter into list + // get indices, physical location and + // boundary values of dofs on this + // face + face->get_dof_indices (face_dofs); + fe_values.reinit(cell, face_no); + const std::vector > &dof_locations = fe_values.get_quadrature_points (); - for (unsigned int i=0; i void VectorTools::project_boundary_values (const DoFHandler &dof, - const typename std::map*> &boundary_functions, + const std::map*> &boundary_functions, const Quadrature &q, std::map &boundary_values) { @@ -858,19 +813,19 @@ VectorTools::integrate_difference (const DoFHandler &dof, if ((norm==H1_seminorm) || (norm==H1_norm)) update_flags = UpdateFlags (update_flags | update_gradients); - FEValues fe_values(fe, q, update_flags); + FEValues fe_values(mapping, fe, q, update_flags); std::vector< Vector > function_values (n_q_points, - Vector(n_components)); + Vector(n_components)); std::vector > > function_grads (n_q_points, - std::vector >(n_components)); + std::vector >(n_components)); std::vector weight_values (n_q_points); std::vector > weight_vectors (n_q_points, n_components); std::vector > psi_values (n_q_points, - Vector(n_components)); + Vector(n_components)); std::vector > > psi_grads (n_q_points, - std::vector >(n_components)); + std::vector >(n_components)); std::vector psi_scalar (n_q_points); // tmp vector when we use the // Function functions for @@ -913,15 +868,15 @@ VectorTools::integrate_difference (const DoFHandler &dof, // the function really has only // one component if (fe_is_system) - exact_solution.vector_value_list (fe_values.get_quadrature_points(), - psi_values); + exact_solution.vector_value_list (fe_values.get_quadrature_points(), + psi_values); else { exact_solution.value_list (fe_values.get_quadrature_points(), tmp_values); for (unsigned int i=0; i &dof, psi_grads[q][k] -= function_grads[q][k]; // take square of integrand - fill_n (psi_scalar.begin(), n_q_points, 0.0); + std::fill_n (psi_scalar.begin(), n_q_points, 0.0); for (unsigned int k=0; k &dof, // add seminorm to L_2 norm or // to zero - diff += inner_product (psi_scalar.begin(), psi_scalar.end(), - fe_values.get_JxW_values().begin(), - 0.0); + diff += std::inner_product (psi_scalar.begin(), psi_scalar.end(), + fe_values.get_JxW_values().begin(), + 0.0); diff = sqrt(diff); break; @@ -1165,20 +1120,18 @@ VectorTools::compute_mean_value (const DoFHandler &dof, Assert (component < dof.get_fe().n_components(), ExcIndexRange(component, 0, dof.get_fe().n_components())); - FEValues fe(dof.get_fe(), quadrature, + FEValues fe(mapping, dof.get_fe(), quadrature, UpdateFlags(update_JxW_values | update_values)); - DoFHandler::active_cell_iterator c; + typename DoFHandler::active_cell_iterator c; std::vector > values(quadrature.n_quadrature_points, - Vector (dof.get_fe().n_components())); + Vector (dof.get_fe().n_components())); double mean = 0.; double area = 0.; // Compute mean value - for (c = dof.begin_active(); - c != dof.end(); - ++c) + for (c = dof.begin_active(); c != dof.end(); ++c) { fe.reinit (c); fe.get_function_values(v, values); @@ -1186,10 +1139,13 @@ VectorTools::compute_mean_value (const DoFHandler &dof, { mean += fe.JxW(k) * values[k](component); area += fe.JxW(k); - } - } + }; + }; + return (mean/area); -} +}; + + @@ -1231,8 +1187,6 @@ double VectorTools::compute_mean_value (const DoFHandler &dof const unsigned int component); - - // the following two functions are not derived from a template in 1d // and thus need no explicit instantiation #if deal_II_dimension > 1 diff --git a/deal.II/examples/step-2/step-2.cc b/deal.II/examples/step-2/step-2.cc index 9cec48bc60..02f9183595 100644 --- a/deal.II/examples/step-2/step-2.cc +++ b/deal.II/examples/step-2/step-2.cc @@ -15,6 +15,7 @@ // ("DoF"s) to vertices, lines, and // cells. #include + // The following include contains the // description of the bilinear finite // element, including the facts that @@ -30,7 +31,7 @@ // elements, but not only for two // space dimensions, but also for one // and three dimensions. -#include +#include // In the following file, several // tools for manipulating degrees of // freedom can be found: @@ -144,21 +145,33 @@ void distribute_dofs (DoFHandler<2> &dof_handler) // object which describes how many // degrees of freedom are to be // associated to each of these - // objects. For (bi-, tri-)linear - // finite elements, this is done - // using the FEQ1 class, which + // objects. This is done using the + // ``FE_Q'' class. Giving a + // constructor argument one // states that one degree of // freedom is to be assigned to // each vertex, while there are // none on lines and inside the - // quadrilateral. We first need to - // create an object of this class - // and use it to distribute the - // degrees of freedom. Note that - // the DoFHandler object will store - // a reference to this object, so - // we need to make it static as - // well, in order to prevent its + // quadrilateral. In fact, the + // argument denotes the polynomial + // degree, in this case, we get + // bilinear finite elements in two + // space dimensions; a value of, + // say, three would give us + // bi-cubic ones. In general, ``FE_Q'' + // denotes the family of continuous + // elements with complete polynomials + // (i.e. tensor-product polynomials) + // up to the specified order + // + // We first need to create an + // object of this class and use it + // to distribute the degrees of + // freedom. Note that the + // DoFHandler object will store a + // reference to this object, so we + // need to make it static as well, + // in order to prevent its // preemptive // destruction. (However, the // library would warn us about this @@ -166,7 +179,7 @@ void distribute_dofs (DoFHandler<2> &dof_handler) // occured. You can check this, if // you want, by removing the // 'static' declaration.) - static const FEQ1<2> finite_element; + static const FE_Q<2> finite_element(1); dof_handler.distribute_dofs (finite_element); // Now we have associated a number diff --git a/deal.II/examples/step-3/step-3.cc b/deal.II/examples/step-3/step-3.cc index 9c6fdc2d4b..f1f7379396 100644 --- a/deal.II/examples/step-3/step-3.cc +++ b/deal.II/examples/step-3/step-3.cc @@ -23,7 +23,7 @@ // In this file are the finite // element descriptions. -#include +#include // And this file is needed for the // creation of sparsity patterns of @@ -118,11 +118,13 @@ class LaplaceProblem // variables. There are variables // describing the triangulation // and the numbering of the - // degrees of freedom... - // (FEQ1<2> is the space with - // shape functions {1,x,y,xy}) + // degrees of freedom (we will + // specify the exact polynomial + // degree of the finite element + // in the constructor of this + // class)... Triangulation<2> triangulation; - FEQ1<2> fe; + FE_Q<2> fe; DoFHandler<2> dof_handler; // ...variables for the sparsity @@ -142,14 +144,20 @@ class LaplaceProblem // Here comes the constructor. It - // does not much more than associate - // the dof_handler variable to the + // does not much more than first to + // specify that we want bi-linear + // elements (denoted by the parameter + // to the finite element object, + // which specifies the polynomial + // degree), and to associate the + // dof_handler variable to the // triangulation we use. All the // other member variables of the // LaplaceProblem class have a // default constructor which does all // we want. LaplaceProblem::LaplaceProblem () : + fe (1), dof_handler (triangulation) {}; diff --git a/deal.II/examples/step-4/step-4.cc b/deal.II/examples/step-4/step-4.cc index 8566ce580d..754ddc6301 100644 --- a/deal.II/examples/step-4/step-4.cc +++ b/deal.II/examples/step-4/step-4.cc @@ -12,7 +12,7 @@ #include #include #include -#include +#include #include #include #include @@ -72,7 +72,7 @@ class LaplaceProblem void output_results () const; Triangulation triangulation; - FEQ1 fe; + FE_Q fe; DoFHandler dof_handler; SparsityPattern sparsity_pattern; @@ -199,12 +199,15 @@ double BoundaryValues::value (const Point &p, // This is the constructor of the - // LaplaceProblem class. It - // associates the DoFHandler to the + // LaplaceProblem class. It specifies + // the desired polynomial degree of + // the finite elements and associates + // the DoFHandler to the // triangulation just as in the // previous example. template LaplaceProblem::LaplaceProblem () : + fe (1), dof_handler (triangulation) {}; @@ -414,8 +417,7 @@ void LaplaceProblem::assemble_system () right_hand_side.value (fe_values.quadrature_point (q_point)) * fe_values.JxW (q_point)); }; - - + // The transfer into the global // matrix and right hand side // is done exactly as before, diff --git a/deal.II/examples/step-5/step-5.cc b/deal.II/examples/step-5/step-5.cc index 5477e710f3..8a6fcb63f7 100644 --- a/deal.II/examples/step-5/step-5.cc +++ b/deal.II/examples/step-5/step-5.cc @@ -19,7 +19,7 @@ #include #include #include -#include +#include #include #include #include @@ -69,7 +69,7 @@ class LaplaceProblem void output_results (const unsigned int cycle) const; Triangulation triangulation; - FEQ1 fe; + FE_Q fe; DoFHandler dof_handler; SparsityPattern sparsity_pattern; @@ -257,6 +257,7 @@ void Coefficient::value_list (const std::vector > &points, // This function is as before. template LaplaceProblem::LaplaceProblem () : + fe (1), dof_handler (triangulation) {}; diff --git a/deal.II/examples/step-6/step-6.cc b/deal.II/examples/step-6/step-6.cc index f617d9843e..3682e23520 100644 --- a/deal.II/examples/step-6/step-6.cc +++ b/deal.II/examples/step-6/step-6.cc @@ -27,14 +27,11 @@ #include #include // From the following include file we - // will import the declaration of the - // quadratic finite element class, - // which in analogy to ``FEQ1'' for - // the linear element is called - // ``FEQ2''. The Lagrange elements of - // polynomial degrees one through four - // are all declared in this file. -#include + // will import the declaration of + // H1-conforming finite element shape + // functions. This family of + // finite elements is called ``FE_Q''. +#include // We will not read the grid from a // file as in the previous example, // but generate it using a function @@ -114,13 +111,14 @@ class LaplaceProblem void output_results (const unsigned int cycle) const; Triangulation triangulation; - DoFHandler dof_handler; - // In order to use the quadratic - // element, we only have to - // replace the declaration of the - // ``fe'' variable like this: - FEQ2 fe; + // We need a finite element + // again. This time, we will want + // to use quadratic polynomials + // (but this is only specified in + // the constructor): + FE_Q fe; + DoFHandler dof_handler; // This is the new variable in // the main class. We need an @@ -187,8 +185,17 @@ void Coefficient::value_list (const std::vector > &points, }; + // This is mostly the same as before, + // but this time we want to use the + // quadratic element. To do so, we + // only have to replace the + // constructor argument (which was + // ``1'' in all previous examples) by + // the desired polynomial degree + // (here ``2''): template LaplaceProblem::LaplaceProblem () : + fe (2), dof_handler (triangulation) {}; -- 2.39.5