// $Id$
// Version: $Name$
//
-// Copyright (C) 1998, 1999, 2000, 2001 by the deal.II authors
+// Copyright (C) 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
* Class for internal data of finite
* element and mapping objects.
*/
+//TODO: can we make the following class protected?
class InternalDataBase: public Subscriptor
{
public:
* Virtual destructor.
*/
virtual ~Mapping ();
-
+
+//TODO: why make the following functions public? they are only helpful for fevalues and maybe the finite elements?
/**
* Prepare internal data
* structures and fill in values
virtual void transform_contravariant (std::vector<Tensor<1,dim> > &dst,
const std::vector<Tensor<1,dim> > &src,
const InternalDataBase& internal,
- const unsigned int src_offset) const = 0;
+ const unsigned int src_offset) const = 0;
/**
* Tranform a field of covariant vectors.
// $Id$
// Version: $Name$
//
-// Copyright (C) 1998, 1999, 2000, 2001 by the deal.II authors
+// Copyright (C) 2001 by the deal.II authors
//
// This file is subject to QPL and may not be distributed
// without copyright and license information. Please refer
// $Id$
// Version: $Name$
//
-// Copyright (C) 1998, 1999, 2000, 2001 by the deal.II authors
+// Copyright (C) 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
*/
~MappingQ ();
+//TODO: why make the following functions public? they are only helpful for fevalues and maybe the finite elements?
/**
* Implementation of the interface in
* @ref{Mapping}.
/**
* Constructor.
*/
- InternalData(unsigned int n_shape_functions);
+ InternalData (const unsigned int n_shape_functions);
/**
* Unit normal vectors. Used
/**
* Takes a
- * @p{set_laplace_on_hex(quad)_vector}
+ * @p{laplace_on_hex(quad)_vector}
* and applies it to the vector
* @p{a} to compute the inner
- * support points. They are
- * appended to the vector @p{a}.
+ * support points as a linear
+ * combination of the exterior
+ * points.
+ *
+ * The vector @p{a} initially
+ * containts the locations of the
+ * @p{n_outer} points, the
+ * @p{n_inner} computed inner
+ * points are appended.
*/
void apply_laplace_vector(const std::vector<std::vector<double> > &lvs,
std::vector<Point<dim> > &a) const;
* solution of a Laplace equation
* with the position of the outer
* support points as boundary
- * values. The outer support
- * points are all support points
- * except of the inner ones.
+ * values, in order to make the
+ * transformation as smooth as
+ * possible.
*/
void compute_support_points_laplace(
const typename Triangulation<dim>::cell_iterator &cell,
std::vector<Point<dim> > &a) const;
/**
- * For @p{dim=2,3}. Adds (appends) the
- * support points of all lines to
- * the vector a.
+ * For @p{dim=2,3}. Append
+ * (appends) the support points
+ * of all shape functions located
+ * on bounding lines to the
+ * vector @p{a}. Points located
+ * on the line but on vertices
+ * are not included.
*
* Needed by the
* @p{compute_support_points_simple(laplace)}
* functions. For @p{dim=1} this
* function is empty.
+ *
+ * This function is made virtual
+ * in order to allow derived
+ * classes to choose shape
+ * function support points
+ * differently than the present
+ * class, which chooses the
+ * points as interpolation points
+ * on the boundary.
*/
- void add_line_support_points (const Triangulation<dim>::cell_iterator &cell,
- std::vector<Point<dim> > &a) const;
+ virtual void
+ add_line_support_points (const Triangulation<dim>::cell_iterator &cell,
+ std::vector<Point<dim> > &a) const;
/**
- * For @p{dim=3}. Adds (appends) the
- * support points of all faces (quads in 3d) to
- * the vector a.
+ * For @p{dim=3}. Append the
+ * support points of all shape
+ * functions located on bounding
+ * faces (quads in 3d) to the
+ * vector @p{a}. Points located
+ * on the line but on vertices
+ * are not included.
*
* Needed by the
* @p{compute_support_points_laplace}
- * function. For @p{dim=1} and 2 this
- * function is empty.
+ * function. For @p{dim=1} and 2
+ * this function is empty.
+ *
+ * This function is made virtual
+ * in order to allow derived
+ * classes to choose shape
+ * function support points
+ * differently than the present
+ * class, which chooses the
+ * points as interpolation points
+ * on the boundary.
*/
//TODO: rename function to add_quad_support_points, to unify notation
- void add_face_support_points(const typename Triangulation<dim>::cell_iterator &cell,
- std::vector<Point<dim> > &a) const;
+ virtual void
+ add_face_support_points(const typename Triangulation<dim>::cell_iterator &cell,
+ std::vector<Point<dim> > &a) const;
/**
* For @p{dim=2} and 3. Simple
* Needed by the
* @p{compute_support_points_simple}
*/
+//TODO: remove this function altogether?
void fill_quad_support_points_simple (const Triangulation<dim>::cell_iterator &cell,
std::vector<Point<dim> > &a) const;
* Number of the Qp tensor
* product shape functions.
*/
- unsigned int n_shape_functions;
+ const unsigned int n_shape_functions;
/**
* Mapping from lexicographic to
// $Id$
// Version: $Name$
//
-// Copyright (C) 1998, 1999, 2000, 2001 by the deal.II authors
+// Copyright (C) 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
class MappingQ1 : public Mapping<dim>
{
public:
+//TODO: why make the following functions public? they are only helpful for fevalues and maybe the finite elements?
/**
* Implementation of the interface in
* @ref{Mapping}.
{
public:
/**
- * Constructor.
+ * Constructor. Pass the
+ * number of shape functions.
*/
- InternalData(unsigned int n_shape_functions);
+ InternalData(const unsigned int n_shape_functions);
/**
* Shape function at quadrature
* vertices must be reordered
* to obtain transformation.
*/
- double shape (unsigned int qpoint,
- unsigned int shape_nr) const;
+ double shape (const unsigned int qpoint,
+ const unsigned int shape_nr) const;
/**
* Shape function at quadrature
* point. See above.
*/
- double &shape (unsigned int qpoint,
- unsigned int shape_nr);
+ double &shape (const unsigned int qpoint,
+ const unsigned int shape_nr);
/**
* Gradient of shape function
* in quadrature point. See
* above.
*/
- Tensor<1,dim> derivative (unsigned int qpoint,
- unsigned int shape_nr) const;
+ Tensor<1,dim> derivative (const unsigned int qpoint,
+ const unsigned int shape_nr) const;
/**
* Gradient of shape function
* in quadrature point. See
* above.
*/
- Tensor<1,dim> &derivative (unsigned int qpoint,
- unsigned int shape_nr);
+ Tensor<1,dim> &derivative (const unsigned int qpoint,
+ const unsigned int shape_nr);
/**
* Values of shape
*/
std::vector<std::vector<Tensor<1,dim> > > aux;
- /**
- * Number of shape functions.
- */
- unsigned int n_shape_functions;
-
/**
* Stores the support points of
* the mapping shape functions on
* @p{mapping_support_points} are
* stored.
*/
- DoFHandler<dim>::cell_iterator cell_of_current_support_points;
+ typename DoFHandler<dim>::cell_iterator cell_of_current_support_points;
/**
* Default value of this flag
* @p{false}.
*/
bool is_mapping_q1_data;
+
+ /**
+ * Number of shape
+ * functions. If this is a Q1
+ * mapping, then it is simply
+ * the number of vertices per
+ * cell. However, since also
+ * derived classes use this
+ * class (e.g. the
+ * @ref{Mapping_Q} class),
+ * the number of shape
+ * functions may also be
+ * different.
+ */
+ unsigned int n_shape_functions;
};
/**
std::vector<Point<dim> > &a) const;
/**
- *Number of shape functions
+ * Number of shape functions. Is
+ * simply the number of vertices
+ * per cell for the Q1 mapping.
*/
- static const unsigned int n_shape_functions = 1 << dim;
+ static const unsigned int n_shape_functions = GeometryInfo<dim>::vertices_per_cell;
};
n_outer((dim==2) ? 4+4*(degree-1)
:8+12*(degree-1)+6*(degree-1)*(degree-1)),
tensor_pols(0),
- n_shape_functions(0),
+ n_shape_functions(power(degree+1,dim)),
renumber(0),
//TODO: why have two ways to compute? if they both work, choose one and remove the other
alternative_normals_computation(false),
v.push_back(LagrangeEquidistant(degree,i));
tensor_pols = new TensorProductPolynomials<dim> (v);
- n_shape_functions=tensor_pols->n_tensor_product_polynomials();
+ Assert (n_shape_functions==tensor_pols->n_tensor_product_polynomials(),
+ ExcInternalError());
Assert(n_inner+n_outer==n_shape_functions, ExcInternalError());
// build the renumbering of the
}
+
template <int dim>
UpdateFlags
MappingQ<dim>::update_each (const UpdateFlags in) const
template <int dim>
void
-MappingQ<dim>::compute_face_data (UpdateFlags update_flags,
- const Quadrature<dim>& q,
- const unsigned int n_original_q_points,
- MappingQ1<dim>::InternalData& mapping_q1_data) const
+MappingQ<dim>::compute_face_data (UpdateFlags update_flags,
+ const Quadrature<dim> &q,
+ const unsigned int n_original_q_points,
+ MappingQ1<dim>::InternalData &mapping_q1_data) const
{
// convert data object to internal
// data for this class. fails with
// possible
InternalData &data = dynamic_cast<InternalData&> (mapping_data);
- data.use_mapping_q1_on_current_cell=!(use_mapping_q_on_all_cells
- || cell->has_boundary_lines());
-
+ // check whether this cell needs
+ // the full mapping or can be
+ // treated by a reduced Q1 mapping,
+ // e.g. if the cell is in the
+ // interior of the domain
+ data.use_mapping_q1_on_current_cell = !(use_mapping_q_on_all_cells
+ || cell->has_boundary_lines());
+
+ // depending on this result, use
+ // this or the other data object
+ // for the mapping
if (data.use_mapping_q1_on_current_cell)
MappingQ1<dim>::fill_fe_values(cell, q, data.mapping_q1_data,
quadrature_points, JxW_values);
}
+
template <int dim>
void
MappingQ<dim>::fill_fe_face_values (const typename DoFHandler<dim>::cell_iterator &cell,
// possible
InternalData &data = dynamic_cast<InternalData&> (mapping_data);
+ // check whether this cell needs
+ // the full mapping or can be
+ // treated by a reduced Q1 mapping,
+ // e.g. if the cell is in the
+ // interior of the domain
//TODO: shouldn't we ask whether the face is at the boundary, rather than the cell?
data.use_mapping_q1_on_current_cell=!(use_mapping_q_on_all_cells
|| cell->has_boundary_lines());
const unsigned int npts=q.n_quadrature_points;
const unsigned int offset=face_no*npts;
+ // depending on this result, use
+ // this or the other data object
+ // for the mapping
if (data.use_mapping_q1_on_current_cell)
MappingQ1<dim>::compute_fill_face (cell, face_no, false,
npts, offset, q.get_weights(),
// possible
InternalData &data = dynamic_cast<InternalData&> (mapping_data);
+ // check whether this cell needs
+ // the full mapping or can be
+ // treated by a reduced Q1 mapping,
+ // e.g. if the cell is in the
+ // interior of the domain
//TODO: shouldn't we ask whether the face is at the boundary, rather than the cell?
data.use_mapping_q1_on_current_cell=!(use_mapping_q_on_all_cells
|| cell->has_boundary_lines());
const unsigned int offset=
(face_no*GeometryInfo<dim>::subfaces_per_face + sub_no)*npts;
+ // depending on this result, use
+ // this or the other data object
+ // for the mapping
if (data.use_mapping_q1_on_current_cell)
MappingQ1<dim>::compute_fill_face (cell, face_no, true,
npts, offset, q.get_weights(),
const unsigned int n_outer_apply=lvs[0].size();
Assert(a.size()==n_outer_apply, ExcDimensionMismatch(a.size(), n_outer_apply));
+ // compute each inner point as
+ // linear combination of the outer
+ // points. the weights are given by
+ // the lvs entries, the outer
+ // points are the first (existing)
+ // elements of a
for (unsigned int unit_point=0; unit_point<n_inner_apply; ++unit_point)
{
- const std::vector<double> &lv=lvs[unit_point];
- Assert(lv.size()==n_outer_apply, ExcInternalError());
+ Assert(lvs[unit_point].size()==n_outer_apply, ExcInternalError());
Point<dim> p;
for (unsigned int k=0; k<n_outer_apply; ++k)
- p+=lv[k]*a[k];
+ p+=lvs[unit_point][k]*a[k];
a.push_back(p);
}
const typename Triangulation<dim>::cell_iterator &cell,
std::vector<Point<dim> > &a) const
{
+ // if this is a cell for which we
+ // want to compute the full
+ // mapping, then get them from the
+ // following function
if (use_mapping_q_on_all_cells || cell->has_boundary_lines())
compute_support_points_laplace(cell, a);
// compute_support_points_simple(cell, a);
+//TODO: can we delete the previous line?
+
else
+ // otherwise: use a Q1 mapping
+ // for which the mapping shape
+ // function support points are
+ // simply the vertices of the
+ // cell
{
a.resize(GeometryInfo<dim>::vertices_per_cell);
MappingQ<dim>::compute_support_points_laplace(const typename Triangulation<dim>::cell_iterator &cell,
std::vector<Point<dim> > &a) const
{
+ // in any case, we need the
+ // vertices first
a.resize(GeometryInfo<dim>::vertices_per_cell);
- // the vertices first
for (unsigned int i=0; i<GeometryInfo<dim>::vertices_per_cell; ++i)
a[i] = cell->vertex(i);
if (degree>1)
- {
- if (dim==1)
- {
- Assert(false, ExcNotImplemented());
- }
- else
- {
- // then the points on lines
- // (for dim=2,3)
- add_line_support_points (cell, a);
-
- if (dim==2)
- apply_laplace_vector(laplace_on_quad_vector,a);
- else if (dim==3)
- {
- add_face_support_points(cell, a);
+ switch (dim)
+ {
+ case 2:
+ // in 2d, add the
+ // points on the four
+ // bounding lines to
+ // the exterior (outer)
+ // points
+ add_line_support_points (cell, a);
+ apply_laplace_vector (laplace_on_quad_vector,a);
+ break;
+
+ case 3:
+ // in 3d also add the
+ // points located on
+ // the boundary faces
+ add_line_support_points (cell, a);
+ add_face_support_points (cell, a);
+ apply_laplace_vector (laplace_on_hex_vector, a);
+ break;
- apply_laplace_vector(laplace_on_hex_vector, a);
- }
- }
- }
+ default 1:
+ Assert(false, ExcNotImplemented());
+ break;
+ };
}
+
+//TODO: remove the following function altogether
template <int dim>
void
MappingQ<dim>::compute_support_points_simple(const typename Triangulation<dim>::cell_iterator &cell,
MappingQ<dim>::add_line_support_points (const Triangulation<dim>::cell_iterator &cell,
std::vector<Point<dim> > &a) const
{
- const Boundary<dim> *boundary = 0;
-
- std::vector<Point<dim> > line_points;
- if (degree>2)
- line_points.resize(degree-1);
-
- // loop over each of the lines, and
- // if it is at the boundary, then
- // first get the boundary
- // description and second compute
- // the points on it
- for (unsigned int line_no=0; line_no<GeometryInfo<dim>::lines_per_cell; ++line_no)
+ // if we only need the midpoint,
+ // then ask for it.
+ if (degree==2)
{
- const typename Triangulation<dim>::line_iterator line = cell->line(line_no);
- if (line->at_boundary())
- boundary=&line->get_triangulation().get_boundary(line->boundary_indicator());
- else
- boundary=&straight_boundary;
-
- // if we only need the
- // midpoint, then ask for
- // it. otherwise call the more
- // complicated functions
- if (degree==2)
+ for (unsigned int line_no=0; line_no<GeometryInfo<dim>::lines_per_cell; ++line_no)
a.push_back(boundary->get_new_point_on_line(line));
- else
+ }
+ else
+ // otherwise call the more
+ // complicated functions and ask
+ // for inner points from the
+ // boundary description
+ {
+ std::vector<Point<dim> > line_points (degree-1);
+
+ // loop over each of the lines,
+ // and if it is at the
+ // boundary, then first get the
+ // boundary description and
+ // second compute the points on
+ // it
+ for (unsigned int line_no=0; line_no<GeometryInfo<dim>::lines_per_cell; ++line_no)
{
+ const typename Triangulation<dim>::line_iterator line = cell->line(line_no);
+
+ const Boundary<dim> * const boundary
+ = (line->at_boundary() ?
+ &line->get_triangulation().get_boundary(line->boundary_indicator()) :
+ &straight_boundary);
+
boundary->get_intermediate_points_on_line (line, line_points);
a.insert (a.end(), line_points.begin(), line_points.end());
- }
+ }
}
-}
+};
template<>
void
MappingQ<3>::add_face_support_points(const Triangulation<3>::cell_iterator &cell,
- std::vector<Point<3> > &a) const
+ std::vector<Point<3> > &a) const
{
const unsigned int faces_per_cell = GeometryInfo<3>::faces_per_cell,
vertices_per_face = GeometryInfo<3>::vertices_per_face,
// on it
if (face->at_boundary())
{
+//TODO: move this variable out of the inner loop
std::vector<Point<3> > quad_points ((degree-1)*(degree-1));
face->get_triangulation().get_boundary(face->boundary_indicator())
// intermediate points
// from a straight
// boundary object
+//TODO: move this variable out of the loop
std::vector<Point<3> > quad_points ((degree-1)*(degree-1));
straight_boundary.get_intermediate_points_on_quad (face, quad_points);
+template <int dim>
+const unsigned int MappingQ1<dim>::n_shape_functions;
+
+
+
template<int dim>
-MappingQ1<dim>::InternalData::InternalData (unsigned int n_shape_functions):
- n_shape_functions(n_shape_functions),
- is_mapping_q1_data(true)
+MappingQ1<dim>::InternalData::InternalData (const unsigned int n_shape_functions):
+ is_mapping_q1_data(true),
+ n_shape_functions (n_shape_functions)
{}
template<int dim> inline
double
-MappingQ1<dim>::InternalData::shape (unsigned int qpoint,
- unsigned int shape_nr) const
+MappingQ1<dim>::InternalData::shape (const unsigned int qpoint,
+ const unsigned int shape_nr) const
{
Assert(qpoint*n_shape_functions + shape_nr < shape_values.size(),
ExcIndexRange(qpoint*n_shape_functions + shape_nr, 0, shape_values.size()));
template <int dim>
void
-MappingQ1<dim>::compute_data (const UpdateFlags update_flags,
- const Quadrature<dim>& q,
- const unsigned int n_original_q_points,
- InternalData& data) const
+MappingQ1<dim>::compute_data (const UpdateFlags update_flags,
+ const Quadrature<dim> &q,
+ const unsigned int n_original_q_points,
+ InternalData &data) const
{
const unsigned int npts = q.n_quadrature_points;
const UpdateFlags flags(data.update_flags);
- // cerr << "Data: " << hex << flags << dec << endl;
-
if (flags & update_transformation_values)
data.shape_values.resize(data.n_shape_functions * npts);
template <int dim>
Mapping<dim>::InternalDataBase*
-MappingQ1<dim>::get_face_data (const UpdateFlags update_flags,
- const Quadrature<dim-1>& quadrature) const
+MappingQ1<dim>::get_face_data (const UpdateFlags update_flags,
+ const Quadrature<dim-1> &quadrature) const
{
InternalData* data = new InternalData(n_shape_functions);
QProjector<dim> q (quadrature, false);
typename std::vector<tensor_>::const_iterator vec = src.begin() + src_offset;
typename std::vector<Tensor<2,dim> >::const_iterator tensor = data.covariant.begin();
typename std::vector<tensor_>::iterator result = dst.begin();
- typename std::vector<tensor_>::const_iterator end = dst.end();
+ const typename std::vector<tensor_>::const_iterator end = dst.end();
while (result!=end)
{