#include <string>
template <int dim> class Triangulation;
+template <int dim> class Mapping;
/**
* This class provides a means to output a triangulation to a file in different
* 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);
};
/**
* 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.
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);
};
* listed in the base class
*/
template <int dim>
- 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.
/**
* 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
* 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
*/
template <int dim>
void write_gnuplot (const Triangulation<dim> &tria,
- std::ostream &out);
+ std::ostream &out,
+ const Mapping<dim> *mapping=0);
/**
* Write the triangulation in the
* 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
*/
template <int dim>
void write_eps (const Triangulation<dim> &tria,
- std::ostream &out);
+ std::ostream &out,
+ const Mapping<dim> *mapping=0);
/**
* Write data and grid to @p{out} according
template <int dim>
void write (const Triangulation<dim> &tria,
std::ostream &out,
- const OutputFormat output_format);
+ const OutputFormat output_format,
+ const Mapping<dim> *mapping=0);
/**
* Set the flags to be used for output
//---------------------------- 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
/**
* 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.
* 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.
*/
Vector<double> &rhs_vector,
const Function<dim> *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<dim> &dof,
- SparseMatrix<double> &matrix);
-
/**
* Assemble the mass matrix and a right
* hand side vector along the boundary.
{}
-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)
#include <base/point.h>
+#include <base/quadrature.h>
#include <grid/grid_out.h>
#include <grid/tria.h>
#include <grid/tria_accessor.h>
#include <grid/tria_iterator.h>
+#include <fe/mapping.h>
#include <iomanip>
#include <algorithm>
};
};
+#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<unsigned int>(cell->material_id()) << std::endl
+ << cell->vertex(1)
+ << ' ' << cell->level()
+ << ' ' << static_cast<unsigned int>(cell->material_id()) << std::endl
+ << std::endl;
+ break;
+ }
+ AssertThrow (out, ExcIO());
+}
+
+
+#else
template <int dim>
void GridOut::write_gnuplot (const Triangulation<dim> &tria,
- std::ostream &out)
+ std::ostream &out,
+ const Mapping<dim> *mapping)
{
AssertThrow (out, ExcIO());
+ const unsigned int n_points=
+ gnuplot_flags.n_boundary_face_points;
+
typename Triangulation<dim>::active_cell_iterator cell=tria.begin_active();
const typename Triangulation<dim>::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<dim> *q_projector=0;
+ if (mapping!=0)
+ {
+ typename std::vector<Point<dim-1> > boundary_points(n_points);
+ for (unsigned int i=0; i<n_points; ++i)
+ boundary_points[i](0)= 1.*(i+1)/(n_points+1);
+
+ Quadrature<dim-1> quadrature(boundary_points, std::vector<double> (n_points, 1));
+ q_projector = new QProjector<dim> (quadrature, false);
+ }
+
for (; cell!=endc; ++cell)
{
if (gnuplot_flags.write_cell_numbers)
switch (dim)
{
case 1:
- out << cell->vertex(0)
- << ' ' << cell->level()
- << ' ' << static_cast<unsigned int>(cell->material_id()) << std::endl
- << cell->vertex(1)
- << ' ' << cell->level()
- << ' ' << static_cast<unsigned int>(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<unsigned int>(cell->material_id()) << std::endl
<< ' ' << static_cast<unsigned int>(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<GeometryInfo<dim>::faces_per_cell; ++face_no)
+ {
+ const typename Triangulation<dim>::face_iterator
+ face = cell->face(face_no);
+ if (face->at_boundary())
+ {
+ out << face->vertex(0)
+ << ' ' << cell->level()
+ << ' ' << static_cast<unsigned int>(cell->material_id())
+ << std::endl;
+
+ const unsigned int offset=face_no*n_points;
+ for (unsigned int i=0; i<n_points; ++i)
+ out << (mapping->transform_unit_to_real_cell
+ (cell, q_projector->point(offset+i)))
+ << ' ' << cell->level()
+ << ' ' << static_cast<unsigned int>(cell->material_id())
+ << std::endl;
+
+ out << face->vertex(1)
+ << ' ' << cell->level()
+ << ' ' << static_cast<unsigned int>(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<unsigned int>(cell->material_id())
+ << std::endl
+ << face->vertex(1)
+ << ' ' << cell->level()
+ << ' ' << static_cast<unsigned int>(cell->material_id())
+ << std::endl
+ << std::endl
+ << std::endl;
+ };
+ };
+ };
+
+ break;
+ };
+
case 3:
- // front face
- out << cell->vertex(0)
- << ' ' << cell->level()
- << ' ' << static_cast<unsigned int>(cell->material_id()) << std::endl
- << cell->vertex(1)
- << ' ' << cell->level()
- << ' ' << static_cast<unsigned int>(cell->material_id()) << std::endl
- << cell->vertex(2)
- << ' ' << cell->level()
- << ' ' << static_cast<unsigned int>(cell->material_id()) << std::endl
- << cell->vertex(3)
- << ' ' << cell->level()
- << ' ' << static_cast<unsigned int>(cell->material_id()) << std::endl
- << cell->vertex(0)
- << ' ' << cell->level()
- << ' ' << static_cast<unsigned int>(cell->material_id()) << std::endl
- << std::endl;
- // back face
- out << cell->vertex(4)
- << ' ' << cell->level()
- << ' ' << static_cast<unsigned int>(cell->material_id()) << std::endl
- << cell->vertex(5)
- << ' ' << cell->level()
- << ' ' << static_cast<unsigned int>(cell->material_id()) << std::endl
- << cell->vertex(6)
- << ' ' << cell->level()
- << ' ' << static_cast<unsigned int>(cell->material_id()) << std::endl
- << cell->vertex(7)
- << ' ' << cell->level()
- << ' ' << static_cast<unsigned int>(cell->material_id()) << std::endl
- << cell->vertex(4)
- << ' ' << cell->level()
- << ' ' << static_cast<unsigned int>(cell->material_id()) << std::endl
- << std::endl;
-
- // now for the four connecting lines
- out << cell->vertex(0)
- << ' ' << cell->level()
- << ' ' << static_cast<unsigned int>(cell->material_id()) << std::endl
- << cell->vertex(4)
- << ' ' << cell->level()
- << ' ' << static_cast<unsigned int>(cell->material_id()) << std::endl
- << std::endl;
- out << cell->vertex(1)
- << ' ' << cell->level()
- << ' ' << static_cast<unsigned int>(cell->material_id()) << std::endl
- << cell->vertex(5)
- << ' ' << cell->level()
- << ' ' << static_cast<unsigned int>(cell->material_id()) << std::endl
- << std::endl;
- out << cell->vertex(2)
- << ' ' << cell->level()
- << ' ' << static_cast<unsigned int>(cell->material_id()) << std::endl
- << cell->vertex(6)
- << ' ' << cell->level()
- << ' ' << static_cast<unsigned int>(cell->material_id()) << std::endl
- << std::endl;
- out << cell->vertex(3)
- << ' ' << cell->level()
- << ' ' << static_cast<unsigned int>(cell->material_id()) << std::endl
- << cell->vertex(7)
- << ' ' << cell->level()
- << ' ' << static_cast<unsigned int>(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<unsigned int>(cell->material_id()) << std::endl
+ << cell->vertex(1)
+ << ' ' << cell->level()
+ << ' ' << static_cast<unsigned int>(cell->material_id()) << std::endl
+ << cell->vertex(2)
+ << ' ' << cell->level()
+ << ' ' << static_cast<unsigned int>(cell->material_id()) << std::endl
+ << cell->vertex(3)
+ << ' ' << cell->level()
+ << ' ' << static_cast<unsigned int>(cell->material_id()) << std::endl
+ << cell->vertex(0)
+ << ' ' << cell->level()
+ << ' ' << static_cast<unsigned int>(cell->material_id()) << std::endl
+ << std::endl;
+ // back face
+ out << cell->vertex(4)
+ << ' ' << cell->level()
+ << ' ' << static_cast<unsigned int>(cell->material_id()) << std::endl
+ << cell->vertex(5)
+ << ' ' << cell->level()
+ << ' ' << static_cast<unsigned int>(cell->material_id()) << std::endl
+ << cell->vertex(6)
+ << ' ' << cell->level()
+ << ' ' << static_cast<unsigned int>(cell->material_id()) << std::endl
+ << cell->vertex(7)
+ << ' ' << cell->level()
+ << ' ' << static_cast<unsigned int>(cell->material_id()) << std::endl
+ << cell->vertex(4)
+ << ' ' << cell->level()
+ << ' ' << static_cast<unsigned int>(cell->material_id()) << std::endl
+ << std::endl;
+
+ // now for the four connecting lines
+ out << cell->vertex(0)
+ << ' ' << cell->level()
+ << ' ' << static_cast<unsigned int>(cell->material_id()) << std::endl
+ << cell->vertex(4)
+ << ' ' << cell->level()
+ << ' ' << static_cast<unsigned int>(cell->material_id()) << std::endl
+ << std::endl;
+ out << cell->vertex(1)
+ << ' ' << cell->level()
+ << ' ' << static_cast<unsigned int>(cell->material_id()) << std::endl
+ << cell->vertex(5)
+ << ' ' << cell->level()
+ << ' ' << static_cast<unsigned int>(cell->material_id()) << std::endl
+ << std::endl;
+ out << cell->vertex(2)
+ << ' ' << cell->level()
+ << ' ' << static_cast<unsigned int>(cell->material_id()) << std::endl
+ << cell->vertex(6)
+ << ' ' << cell->level()
+ << ' ' << static_cast<unsigned int>(cell->material_id()) << std::endl
+ << std::endl;
+ out << cell->vertex(3)
+ << ' ' << cell->level()
+ << ' ' << static_cast<unsigned int>(cell->material_id()) << std::endl
+ << cell->vertex(7)
+ << ' ' << cell->level()
+ << ' ' << static_cast<unsigned int>(cell->material_id()) << std::endl
+ << std::endl;
+ break;
+ };
};
};
+
+ if (q_projector != 0)
+ delete q_projector;
+
AssertThrow (out, ExcIO());
};
-
+#endif
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 <int dim>
void GridOut::write_eps (const Triangulation<dim> &tria,
- std::ostream &out)
+ std::ostream &out,
+ const Mapping<dim> *mapping)
{
typedef std::list<LineEntry> LineList;
// 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<EpsFlagsBase&>(eps_flags_2) :
+ (dim==3 ?
+ static_cast<EpsFlagsBase&>(eps_flags_3) :
+ *static_cast<EpsFlagsBase*>(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
switch (dim)
{
+ case 1:
+ {
+ Assert(false, ExcInternalError());
+ break;
+ };
+
case 2:
{
Triangulation<dim>::active_line_iterator line =tria.begin_active_line ();
Triangulation<dim>::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<dim>,
- // but we want a Point<2>.
- // in fact, whenever we're here,
- // the vertex is a Point<dim>,
- // 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<dim>,
+ // but we want a Point<2>.
+ // in fact, whenever we're here,
+ // the vertex is a Point<dim>,
+ // 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<Point<dim-1> > boundary_points (n_points);
+
+ for (unsigned int i=0; i<n_points; ++i)
+ boundary_points[i](0) = 1.*(i+1)/(n_points+1);
+
+ Quadrature<dim-1> quadrature (boundary_points,
+ std::vector<double> (n_points, 1));
+ QProjector<dim> q_projector (quadrature, false);
+
+ // next loop over all
+ // boundary faces and
+ // generate the info from
+ // them
+ typename Triangulation<dim>::active_cell_iterator cell=tria.begin_active ();
+ const typename Triangulation<dim>::active_cell_iterator end=tria.end ();
+ for (; cell!=end; ++cell)
+ for (unsigned int face_no=0; face_no<GeometryInfo<dim>::faces_per_cell; ++face_no)
+ {
+ const typename Triangulation<dim>::face_iterator
+ face = cell->face(face_no);
+
+ if (face->at_boundary())
+ {
+ Point<dim> 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<n_points; ++i)
+ {
+ const Point<dim> 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<dim> 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<dim>::active_line_iterator line =tria.begin_active_line ();
Triangulation<dim>::active_line_iterator endline=tria.end_line ();
// 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<dim> 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.),
};
+
// find out minimum and maximum x and
// y coordinates to compute offsets
// and scaling factors
// lower left corner
<< "0 0 "
// upper right corner
- << static_cast<unsigned int>( (x_max-x_min) * scale )
+ << static_cast<unsigned int>( (x_max-x_min) * scale )+1
<< ' '
- << static_cast<unsigned int>( (y_max-y_min) * scale )
+ << static_cast<unsigned int>( (y_max-y_min) * scale )+1
<< std::endl;
// define some abbreviations to keep
AssertThrow (out, ExcIO());
};
+#endif
template <int dim>
void GridOut::write (const Triangulation<dim> &tria,
std::ostream &out,
- OutputFormat output_format)
+ OutputFormat output_format,
+ const Mapping<dim> *mapping)
{
switch (output_format)
{
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;
};
// explicit instantiations
template void GridOut::write_ucd (const Triangulation<deal_II_dimension> &,
std::ostream &);
+#if deal_II_dimension != 1
template void GridOut::write_gnuplot (const Triangulation<deal_II_dimension> &,
- std::ostream &);
+ std::ostream &,
+ const Mapping<deal_II_dimension> *);
template void GridOut::write_eps (const Triangulation<deal_II_dimension> &,
- std::ostream &);
+ std::ostream &,
+ const Mapping<deal_II_dimension> *);
+#endif
template void GridOut::write (const Triangulation<deal_II_dimension> &,
std::ostream &,
OutputFormat);
#include <lac/full_matrix.h>
#include <lac/vector.h>
#include <lac/sparse_matrix.h>
+#include <base/quadrature.h>
+#include <fe/mapping_q1.h>
+
// if necessary try to work around a bug in the IBM xlC compiler
#ifdef XLC_WORK_AROUND_STD_BUG
#endif
+
+static MappingQ1<deal_II_dimension> mapping;
+
template <int dim>
Assembler<dim>::AssemblerData::AssemblerData (const DoFHandler<dim> &dof,
const bool assemble_matrix,
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)
{
#include <algorithm>
#include <cmath>
-// if necessary try to work around a bug in the IBM xlC compiler
-#ifdef XLC_WORK_AROUND_STD_BUG
-using namespace std;
-#endif
-
template <int dim>
ProblemBase<dim>::ProblemBase () :
template <int dim>
-void ProblemBase<dim>::fill_data (DataOut<dim> &out) const {
+void ProblemBase<dim>::fill_data (DataOut<dim> &out) const
+{
Assert ((tria!=0) && (dof_handler!=0), ExcNoTriaSelected());
out.clear ();
template <int dim>
-std::string ProblemBase<dim>::get_solution_name () const {
+std::string ProblemBase<dim>::get_solution_name () const
+{
return "solution";
};
#include <grid/tria_iterator.h>
#include <fe/fe.h>
#include <fe/fe_values.h>
+#include <fe/mapping_q1.h>
+
+//TODO: Do this more clever
+static MappingQ1<deal_II_dimension> mapping;
#ifdef DEAL_II_USE_MT
#include <base/thread_management.h>
QTrapez<1> q_trapez;
QIterated<dim> patch_points (q_trapez, data.n_subdivisions);
- FEValues<dim> fe_patch_values(dofs->get_fe(),
+ FEValues<dim> fe_patch_values(mapping, dofs->get_fe(),
patch_points,
update_values);
#include <fe/fe.h>
#include <fe/fe_values.h>
+#include <fe/mapping_q1.h>
+
+//TODO: Do this more clever
+static MappingQ1<deal_II_dimension> mapping;
+
#ifdef DEAL_II_USE_MT
#include <base/thread_management.h>
#endif
QTrapez<1> q_trapez;
QIterated<dim-1> patch_points (q_trapez, data.n_subdivisions);
- FEFaceValues<dim> fe_patch_values(dofs->get_fe(),
+ FEFaceValues<dim> fe_patch_values(mapping,
+ dofs->get_fe(),
patch_points,
update_values);
#include <dofs/dof_handler.h>
#include <dofs/dof_accessor.h>
#include <grid/tria_iterator.h>
+#include <fe/mapping_q1.h>
#include <fe/fe.h>
#include <fe/fe_values.h>
+//TODO: Do this more clever
+static MappingQ1<deal_II_dimension> mapping;
+
#ifdef DEAL_II_USE_MT
#include <base/thread_management.h>
#endif
#endif
-
-
template <int dim>
void DataOutRotation<dim>::build_some_patches (Data data)
{
QTrapez<1> q_trapez;
QIterated<dim> patch_points (q_trapez, data.n_subdivisions);
- FEValues<dim> fe_patch_values(dofs->get_fe(),
+ FEValues<dim> fe_patch_values(mapping,
+ dofs->get_fe(),
patch_points,
update_values);
#include <grid/tria_iterator.h>
#include <fe/fe.h>
#include <fe/fe_values.h>
+#include <fe/mapping_q1.h>
+
// if necessary try to work around a bug in the IBM xlC compiler
#endif
+//TODO: Comment?? Use proper mapping!
+static MappingQ1<deal_II_dimension> mapping;
+
template <int dim>
unsigned int
template <int dim>
-void DataOutStack<dim>::declare_data_vector (const std::string &name,
- const VectorType vector_type)
+void DataOutStack<dim>::declare_data_vector (const std::string &name,
+ const VectorType vector_type)
{
std::vector<std::string> names;
names.push_back (name);
template <int dim>
void DataOutStack<dim>::declare_data_vector (const std::vector<std::string> &names,
- const VectorType vector_type)
+ const VectorType vector_type)
{
// make sure this function is
// not called after some parameter
template <int dim>
template <typename number>
void DataOutStack<dim>::add_data_vector (const Vector<number> &vec,
- const std::string &name)
+ const std::string &name)
{
std::vector<std::string> names;
names.push_back (name);
{
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(),
// cell to these points
QTrapez<1> q_trapez;
QIterated<dim> patch_points (q_trapez, n_subdivisions);
- FEValues<dim> fe_patch_values (dof_handler->get_fe(),
+ FEValues<dim> 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<double> patch_values (n_q_points);
std::vector<Vector<double> > patch_values_system (n_q_points,
- Vector<double>(n_components));
+ Vector<double>(n_components));
// add the required number of patches
DataOutBase::Patch<dim+1> default_patch;
// explicit instantiations
template class DataOutStack<deal_II_dimension>;
template void DataOutStack<deal_II_dimension>::add_data_vector (const Vector<double> &vec,
- const std::string &name);
+ const std::string &name);
template void DataOutStack<deal_II_dimension>::add_data_vector (const Vector<float> &vec,
- const std::string &name);
+ const std::string &name);
-//---------------------------- derivative_approximation.cc ---------------------------
+//---------------------------- gradient_estimator.cc ---------------------------
// $Id$
// Version: $Name$
//
#include <dofs/dof_handler.h>
#include <fe/fe.h>
#include <fe/fe_values.h>
+#include <fe/mapping_q1.h>
#include <numerics/derivative_approximation.h>
#include <cmath>
+//TODO: Comment? Use proper mapping!
+static MappingQ1<deal_II_dimension> mapping;
+
template <typename T>
static T sqr (const T t)
Vector<float> &derivative_norm)
{
QMidpoint<dim> midpoint_rule;
- FEValues<dim> fe_midpoint_value (dof_handler.get_fe(),
+ FEValues<dim> fe_midpoint_value (mapping, dof_handler.get_fe(),
midpoint_rule,
UpdateFlags(DerivativeDescription::update_flags |
update_q_points));
//---------------------------- error_estimator.cc ---------------------------
-#include <fe/fe.h>
-#include <fe/fe_values.h>
-#include <fe/fe_update_flags.h>
+#include <base/timer.h>
#include <base/thread_management.h>
#include <base/quadrature.h>
#include <base/quadrature_lib.h>
-#include <numerics/error_estimator.h>
-#include <dofs/dof_handler.h>
+#include <lac/vector.h>
#include <grid/tria_iterator.h>
-#include <dofs/dof_accessor.h>
#include <grid/geometry_info.h>
-#include <lac/vector.h>
+#include <dofs/dof_handler.h>
+#include <dofs/dof_accessor.h>
+#include <fe/fe.h>
+#include <fe/fe_values.h>
+#include <fe/fe_update_flags.h>
+#include <fe/mapping_q1.h>
+#include <numerics/error_estimator.h>
#include <numeric>
#include <algorithm>
#include <cmath>
#include <vector>
-#include <base/timer.h>
-
// if necessary try to work around a bug in the IBM xlC compiler
#ifdef XLC_WORK_AROUND_STD_BUG
#endif
+
+//TODO: Comment?? Proper Mapping
+static MappingQ1<deal_II_dimension> mapping;
+
+
+
static
inline
double sqr (const double x)
};
-
-
#if deal_II_dimension == 1
template <>
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
const FunctionMap &neumann_bc,
const Vector<double> &solution,
Vector<float> &error,
- const std::vector<bool> &component_mask,
+ const std::vector<bool> &component_mask,
const Function<dim> *coefficients,
unsigned int n_threads)
{
// if no mask given: treat all components
std::vector<bool> component_mask ((component_mask_.size() == 0) ?
- std::vector<bool>(n_components, true) :
- component_mask_);
+ std::vector<bool>(n_components, true) :
+ component_mask_);
Assert (component_mask.size() == n_components, ExcInvalidComponentMask());
Assert (count(component_mask.begin(), component_mask.end(), true) > 0,
ExcInvalidComponentMask());
// 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)
{
// need not compute all values on the
// neighbor cells, so using two objects
// gives us a performance gain).
- FEFaceValues<dim> fe_face_values_cell (data.dof_handler.get_fe(),
+ FEFaceValues<dim> fe_face_values_cell (mapping,
+ data.dof_handler.get_fe(),
data.quadrature,
UpdateFlags(update_gradients |
update_JxW_values |
(data.coefficients != 0)) ?
update_q_points : 0) |
update_normal_vectors));
- FEFaceValues<dim> fe_face_values_neighbor (data.dof_handler.get_fe(),
+ FEFaceValues<dim> fe_face_values_neighbor (mapping,
+ data.dof_handler.get_fe(),
data.quadrature,
update_gradients);
- FESubfaceValues<dim> fe_subface_values (data.dof_handler.get_fe(),
+ FESubfaceValues<dim> fe_subface_values (mapping,
+ data.dof_handler.get_fe(),
data.quadrature,
update_gradients);
solutions,
((component_mask.size() == 0) ?
std::vector<bool>(dof_handler.get_fe().n_components(),
- true) :
+ true) :
component_mask),
coefficients,
n_threads,
data.face_integrals.end(),
ExcInternalError());
Assert (data.face_integrals[face->child(subface_no)][0] >= 0,
- ExcInternalError());
-
+ ExcInternalError());
+
for (unsigned int n=0; n<n_solution_vectors; ++n)
sum[n] += data.face_integrals[face->child(subface_no)][n];
};
#include <lac/block_vector.h>
#include <lac/sparse_matrix.h>
#include <lac/block_sparse_matrix.h>
+#include <fe/mapping_q1.h>
#include <algorithm>
#include <set>
#include <cmath>
+
// 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<deal_II_dimension> mapping;
+
+
+
+//TODO: re-create the create_mass_matrix function with 2 args
+
template <int dim>
void MatrixCreator<dim>::create_mass_matrix (const DoFHandler<dim> &dof,
const Quadrature<dim> &q,
-template <int dim>
-void MatrixCreator<dim>::create_mass_matrix (const DoFHandler<dim> &dof,
- SparseMatrix<double> &matrix)
-{
- const FiniteElement<dim> &fe = dof.get_fe();
-
- const unsigned int dofs_per_cell = fe.dofs_per_cell;
-
- FullMatrix<double> local_mass_matrix (dofs_per_cell, dofs_per_cell);
- std::vector<unsigned int> dofs_on_this_cell (dofs_per_cell);
-
- DoFHandler<dim>::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<dofs_per_cell; ++i)
- for (unsigned int j=0; j<dofs_per_cell; ++j)
- matrix.add (dofs_on_this_cell[i], dofs_on_this_cell[j],
- local_mass_matrix(i,j));
- };
-};
-
-
#if deal_II_dimension == 1
template <>
-void MatrixCreator<1>::create_boundary_mass_matrix (const DoFHandler<1> &,
- const Quadrature<0> &,
- SparseMatrix<double> &,
- const FunctionMap &,
- Vector<double> &,
- std::vector<unsigned int> &,
- const Function<1> *)
+void MatrixCreator<1>::create_boundary_mass_matrix (const DoFHandler<1> &,
+ const Quadrature<0> &,
+ SparseMatrix<double> &,
+ const FunctionMap &,
+ Vector<double> &,
+ std::vector<unsigned int> &,
+ const Function<1> *)
{
Assert (false, ExcNotImplemented());
};
template <int dim>
-void MatrixCreator<dim>::create_boundary_mass_matrix (const DoFHandler<dim> &dof,
- const Quadrature<dim-1> &q,
- SparseMatrix<double> &matrix,
- const FunctionMap &rhs,
- Vector<double> &rhs_vector,
- std::vector<unsigned int> &dof_to_boundary_mapping,
- const Function<dim> *a)
+void MatrixCreator<dim>::create_boundary_mass_matrix (const DoFHandler<dim> &dof,
+ const Quadrature<dim-1> &q,
+ SparseMatrix<double> &matrix,
+ const FunctionMap &rhs,
+ Vector<double> &rhs_vector,
+ std::vector<unsigned int> &dof_to_boundary_mapping,
+ const Function<dim> *a)
{
const FiniteElement<dim> &fe = dof.get_fe();
const unsigned int n_components = fe.n_components();
UpdateFlags update_flags = UpdateFlags (update_values |
update_JxW_values |
update_q_points);
- FEFaceValues<dim> fe_values (fe, q, update_flags);
+ FEFaceValues<dim> fe_values (mapping, fe, q, update_flags);
// two variables for the coefficient,
// one for the two cases indicated in
+
+//TODO: recreate this function
/*
template <int dim>
+
+
template <int dim>
void MatrixCreator<dim>::create_laplace_matrix (const DoFHandler<dim> &dof,
const Quadrature<dim> &q,
template <typename number>
void
MatrixTools<dim>::apply_boundary_values (const std::map<unsigned int,double> &boundary_values,
- SparseMatrix<number> &matrix,
- Vector<number> &solution,
- Vector<number> &right_hand_side,
- const bool preserve_symmetry)
+ SparseMatrix<number> &matrix,
+ Vector<number> &solution,
+ Vector<number> &right_hand_side,
+ const bool preserve_symmetry)
{
Assert (matrix.n() == matrix.m(),
ExcDimensionsDontMatch(matrix.n(), matrix.m()));
return;
- std::map<unsigned int,double>::const_iterator dof = boundary_values.begin(),
- endd = boundary_values.end();
+ std::map<unsigned int,double>::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();
BlockSparseMatrix<double> &matrix,
BlockVector<double> &solution,
BlockVector<double> &right_hand_side,
- const bool preserve_symmetry)
+ const bool preserve_symmetry)
{
const unsigned int blocks = matrix.n_block_rows();
return;
- std::map<unsigned int,double>::const_iterator dof = boundary_values.begin(),
- endd = boundary_values.end();
+ std::map<unsigned int,double>::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();
[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
+
+
template <int dim>
LaplaceMatrix<dim>::LaplaceMatrix (const Function<dim> * const rhs,
const Function<dim> * const a) :
#include <lac/precondition.h>
#include <lac/solver_cg.h>
#include <lac/vector_memory.h>
+#include <fe/mapping_q1.h>
#include <numeric>
#include <algorithm>
#endif
+
+//TODO: Comment?? Use proper mapping!
+static MappingQ1<deal_II_dimension> mapping;
+
+
+
static inline double sqr (const double x)
{
return x*x;
+
template <int dim>
void VectorTools::interpolate (const DoFHandler<dim> &dof,
const Function<dim> &function,
// avoided to evaluate
// the vectorfunction multiply at
// the same point on a cell.
- std::vector<Point<dim> > unit_support_points (fe.dofs_per_cell);
+ std::vector<Point<dim> > 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<bool> 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<fe.dofs_per_cell; ++i)
-// {
-// if (unit_support_points[i] != unit_support_points[last_rep_point])
-// {
-// is_representative_point[i] = true;
-// last_rep_point=i;
-// ++n_rep_points;
-// }
-// };
-
-// vector<unsigned int> dofs_on_cell (fe.dofs_per_cell);
-// vector<Point<dim> > support_points (fe.dofs_per_cell);
-
-// vector<Point<dim> > rep_points (n_rep_points);
-// vector<Vector<double> > function_values_at_rep_points (
-// n_rep_points, Vector<double>(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; i<fe.dofs_per_cell; ++i)
-// if (is_representative_point[i])
-// rep_points[j++]=support_points[i];
-// Assert(j == n_rep_points, ExcInternalError());
-
-// // get function values at these points
-// vectorfunction.value_list (rep_points, function_values_at_rep_points);
-
-// // get indices of the dofs on this cell
-// cell->get_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; i<fe.dofs_per_cell; ++i)
-// {
-// if (is_representative_point[i])
-// ++last_rep_point;
-
-// const unsigned int component
-// = fe.system_to_component_index(i).first;
-// vec(dofs_on_cell[i])
-// = function_values_at_rep_points[last_rep_point](component);
-// }
-// }
-
- // The following is more general.
- // It also works if @p{dofs_per_x>1},
- // 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
Assert(dof_to_rep_index_table.size()==fe.dofs_per_cell, ExcInternalError());
std::vector<unsigned int> dofs_on_cell (fe.dofs_per_cell);
- std::vector<Point<dim> > support_points (fe.dofs_per_cell);
+ std::vector<double> dummy_weights (fe.dofs_per_cell);
std::vector<Point<dim> > rep_points (n_rep_points);
std::vector<Vector<double> > function_values_system (n_rep_points,
Vector<double>(fe.n_components()));
+ // Make a quadrature rule from support points
+ // to feed it into FEValues
+ Quadrature<dim> support_quadrature(unit_support_points, dummy_weights);
+
+//TODO: Higher order mapping?
+ // Transformed support points are computed by
+ // FEValues
+ FEValues<dim> 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<Point<dim> >& support_points =
+ fe_values.get_quadrature_points();
// pick out the representative
// support points
+
template <int dim> void
VectorTools::interpolate (const DoFHandler<dim> &dof_1,
const DoFHandler<dim> &dof_2,
SparseMatrix<double> mass_matrix (sparsity);
Vector<double> tmp (mass_matrix.n());
- // try to assemble the mass matrix by exact
- // integration. if this is not supported,
- // then use quadrature
- try
- {
- MatrixCreator<dim>::create_mass_matrix (dof, mass_matrix);
- }
- catch (FiniteElement<dim>::ExcComputationNotUseful)
- {
- MatrixCreator<dim>::create_mass_matrix (dof, quadrature, mass_matrix);
- };
+ MatrixCreator<dim>::create_mass_matrix (dof, quadrature, mass_matrix);
VectorTools::create_right_hand_side (dof, quadrature, function, tmp);
UpdateFlags update_flags = UpdateFlags(update_values |
update_q_points |
update_JxW_values);
- FEValues<dim> fe_values (fe, quadrature, update_flags);
+ 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,
std::vector<unsigned int> dofs (dofs_per_cell);
Vector<double> cell_vector (dofs_per_cell);
- DoFHandler<dim>::active_cell_iterator cell = dof_handler.begin_active(),
- endc = dof_handler.end();
+ typename DoFHandler<dim>::active_cell_iterator cell = dof_handler.begin_active(),
+ endc = dof_handler.end();
if (n_components==1)
{
{
fe_values.reinit(cell);
- const FullMatrix<double> &values = fe_values.get_shape_values ();
- const std::vector<double> &weights = fe_values.get_JxW_values ();
+ const FullMatrix<double> &values = fe_values.get_shape_values ();
+ const std::vector<double> &weights = fe_values.get_JxW_values ();
rhs_function.value_list (fe_values.get_quadrature_points(), rhs_values);
cell_vector.clear();
{
fe_values.reinit(cell);
- const FullMatrix<double> &values = fe_values.get_shape_values ();
- const std::vector<double> &weights = fe_values.get_JxW_values ();
+ const FullMatrix<double> &values = fe_values.get_shape_values ();
+ const std::vector<double> &weights = fe_values.get_JxW_values ();
rhs_function.vector_value_list (fe_values.get_quadrature_points(), rhs_values);
cell_vector.clear();
// the original value or a vector
// of @p{true}s
const std::vector<bool> component_mask ((component_mask_.size() == 0) ?
- std::vector<bool> (fe.n_components(), true) :
- component_mask_);
+ std::vector<bool> (fe.n_components(), true) :
+ component_mask_);
Assert (count(component_mask.begin(), component_mask.end(), true) > 0,
ExcComponentMismatch());
// the original value or a vector
// of @p{true}s
const std::vector<bool> component_mask ((component_mask_.size() == 0) ?
- std::vector<bool> (fe.n_components(), true) :
- component_mask_);
+ std::vector<bool> (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<unsigned int> face_dofs (fe.dofs_per_face, -1);
std::vector<Point<dim> > dof_locations (face_dofs.size(), Point<dim>());
+
// array to store the values of
// the boundary function at the
// boundary points. have to arrays
// respectively
std::vector<double> dof_values_scalar (fe.dofs_per_face);
std::vector<Vector<double> > dof_values_system (fe.dofs_per_face,
- Vector<double>(fe.n_components()));
-
- DoFHandler<dim>::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<double>(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<Point<dim-1> > 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<double> dummy_weights (unit_support_points.size());
+ Quadrature<dim-1> aux_quad (unit_support_points, dummy_weights);
+ FEFaceValues<dim> fe_values (mapping, fe, aux_quad, update_q_points);
+
+ DoFHandler<dim>::active_cell_iterator cell = dof.begin_active(),
+ endc = dof.end();
+ for (; cell!=endc; ++cell)
+ for (unsigned int face_no = 0; face_no < GeometryInfo<dim>::faces_per_cell;
+ ++face_no)
+ {
+ DoFHandler<dim>::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<Point<dim> > &dof_locations = fe_values.get_quadrature_points ();
- for (unsigned int i=0; i<face_dofs.size(); ++i)
- if (component_mask[fe.face_system_to_component_index(i).first])
- boundary_values[face_dofs[i]]
- = dof_values_system[i](fe.face_system_to_component_index(i).first);
+ if (fe_is_system)
+ {
+ boundary_function.vector_value_list (dof_locations, dof_values_system);
+
+ // enter into list
+
+ for (unsigned int i=0; i<face_dofs.size(); ++i)
+ if (component_mask[fe.face_system_to_component_index(i).first])
+ boundary_values[face_dofs[i]]
+ = dof_values_system[i](fe.face_system_to_component_index(i).first);
+ }
+ else
+ // fe has only one component,
+ // so save some computations
+ {
+ // get only the one component that
+ // this function has
+ boundary_function.value_list (dof_locations,
+ dof_values_scalar,
+ 0);
+
+ // enter into list
+
+ for (unsigned int i=0; i<face_dofs.size(); ++i)
+ boundary_values[face_dofs[i]] = dof_values_scalar[i];
+ }
}
- else
- // fe has only one component,
- // so save some computations
- {
- // get only the one component that
- // this function has
- boundary_function.value_list (dof_locations,
- dof_values_scalar,
- 0);
-
- // enter into list
-
- for (unsigned int i=0; i<face_dofs.size(); ++i)
- boundary_values[face_dofs[i]] = dof_values_scalar[i];
- };
- };
+ }
}
-
-
+
+
template <int dim>
void
VectorTools::project_boundary_values (const DoFHandler<dim> &dof,
- const typename std::map<unsigned char,const Function<dim>*> &boundary_functions,
+ const std::map<unsigned char,const Function<dim>*> &boundary_functions,
const Quadrature<dim-1> &q,
std::map<unsigned int,double> &boundary_values)
{
if ((norm==H1_seminorm) || (norm==H1_norm))
update_flags = UpdateFlags (update_flags | update_gradients);
- FEValues<dim> fe_values(fe, q, update_flags);
+ FEValues<dim> fe_values(mapping, fe, q, update_flags);
std::vector< Vector<double> > function_values (n_q_points,
- Vector<double>(n_components));
+ Vector<double>(n_components));
std::vector<std::vector<Tensor<1,dim> > > function_grads (n_q_points,
- std::vector<Tensor<1,dim> >(n_components));
+ std::vector<Tensor<1,dim> >(n_components));
std::vector<double> weight_values (n_q_points);
std::vector<Vector<double> > weight_vectors (n_q_points, n_components);
std::vector<Vector<double> > psi_values (n_q_points,
- Vector<double>(n_components));
+ Vector<double>(n_components));
std::vector<std::vector<Tensor<1,dim> > > psi_grads (n_q_points,
- std::vector<Tensor<1,dim> >(n_components));
+ std::vector<Tensor<1,dim> >(n_components));
std::vector<double> psi_scalar (n_q_points);
// tmp vector when we use the
// Function<dim> functions for
// 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<n_q_points; ++i)
psi_values[i](0) = tmp_values[i];
- }
+ };
// then subtract finite element
// fe_function
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<n_components; ++k)
if (weight != 0)
// 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;
Assert (component < dof.get_fe().n_components(),
ExcIndexRange(component, 0, dof.get_fe().n_components()));
- FEValues<dim> fe(dof.get_fe(), quadrature,
+ FEValues<dim> fe(mapping, dof.get_fe(), quadrature,
UpdateFlags(update_JxW_values
| update_values));
- DoFHandler<dim>::active_cell_iterator c;
+ typename DoFHandler<dim>::active_cell_iterator c;
std::vector<Vector<double> > values(quadrature.n_quadrature_points,
- Vector<double> (dof.get_fe().n_components()));
+ Vector<double> (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);
{
mean += fe.JxW(k) * values[k](component);
area += fe.JxW(k);
- }
- }
+ };
+ };
+
return (mean/area);
-}
+};
+
+
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
// ("DoF"s) to vertices, lines, and
// cells.
#include <dofs/dof_handler.h>
+
// The following include contains the
// description of the bilinear finite
// element, including the facts that
// elements, but not only for two
// space dimensions, but also for one
// and three dimensions.
-#include <fe/fe_lib.lagrange.h>
+#include <fe/fe_q.h>
// In the following file, several
// tools for manipulating degrees of
// freedom can be found:
// 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
// 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
// In this file are the finite
// element descriptions.
-#include <fe/fe_lib.lagrange.h>
+#include <fe/fe_q.h>
// And this file is needed for the
// creation of sparsity patterns of
// 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
// 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)
{};
#include <grid/tria_accessor.h>
#include <grid/tria_iterator.h>
#include <dofs/dof_accessor.h>
-#include <fe/fe_lib.lagrange.h>
+#include <fe/fe_q.h>
#include <dofs/dof_tools.h>
#include <fe/fe_values.h>
#include <base/quadrature_lib.h>
void output_results () const;
Triangulation<dim> triangulation;
- FEQ1<dim> fe;
+ FE_Q<dim> fe;
DoFHandler<dim> dof_handler;
SparsityPattern sparsity_pattern;
// 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 <int dim>
LaplaceProblem<dim>::LaplaceProblem () :
+ fe (1),
dof_handler (triangulation)
{};
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,
#include <grid/tria_iterator.h>
#include <dofs/dof_accessor.h>
#include <dofs/dof_tools.h>
-#include <fe/fe_lib.lagrange.h>
+#include <fe/fe_q.h>
#include <fe/fe_values.h>
#include <numerics/vectors.h>
#include <numerics/matrices.h>
void output_results (const unsigned int cycle) const;
Triangulation<dim> triangulation;
- FEQ1<dim> fe;
+ FE_Q<dim> fe;
DoFHandler<dim> dof_handler;
SparsityPattern sparsity_pattern;
// This function is as before.
template <int dim>
LaplaceProblem<dim>::LaplaceProblem () :
+ fe (1),
dof_handler (triangulation)
{};
#include <numerics/matrices.h>
#include <numerics/data_out.h>
// 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 <fe/fe_lib.lagrange.h>
+ // will import the declaration of
+ // H1-conforming finite element shape
+ // functions. This family of
+ // finite elements is called ``FE_Q''.
+#include <fe/fe_q.h>
// We will not read the grid from a
// file as in the previous example,
// but generate it using a function
void output_results (const unsigned int cycle) const;
Triangulation<dim> triangulation;
- DoFHandler<dim> dof_handler;
- // In order to use the quadratic
- // element, we only have to
- // replace the declaration of the
- // ``fe'' variable like this:
- FEQ2<dim> 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<dim> fe;
+ DoFHandler<dim> dof_handler;
// This is the new variable in
// the main class. We need an
};
+ // 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 <int dim>
LaplaceProblem<dim>::LaplaceProblem () :
+ fe (2),
dof_handler (triangulation)
{};