Mapping<dim,spacedim> *clone () const;
-protected:
-
- /**
- * For <tt>dim=2,3</tt>. Append the support points of all shape functions
- * located on bounding lines to the vector @p a. Points located on the line
- * but not on vertices are not included.
- *
- * Needed by the @p compute_support_points_laplace function . For
- * <tt>dim=1</tt> 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.
- */
- virtual void
- add_line_support_points (const typename Triangulation<dim,spacedim>::cell_iterator &cell,
- std::vector<Point<spacedim> > &a) const;
-
- /**
- * For <tt>dim=3</tt>. Append the support points of all shape functions
- * located on bounding faces (quads in 3d) to the vector @p a. Points
- * located on the quad but not on vertices are not included.
- *
- * Needed by the @p compute_support_points_laplace function. For
- * <tt>dim=1</tt> and <tt>dim=2</tt> 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.
- */
- virtual void
- add_quad_support_points(const typename Triangulation<dim,spacedim>::cell_iterator &cell,
- std::vector<Point<spacedim> > &a) const;
-
-
-private:
- /**
- * Ask the manifold descriptor to return intermediate points on lines or
- * faces. The function needs to return one or multiple points (depending on
- * the number of elements in the output vector @p points that lie inside a
- * line, quad or hex). Whether it is a line, quad or hex doesn't really
- * matter to this function but it can be inferred from the number of input
- * points in the @p surrounding_points vector.
- */
- void get_intermediate_points(const Manifold<dim, spacedim> &manifold,
- const std::vector<Point<spacedim> > &surrounding_points,
- std::vector<Point<spacedim> > &points) const;
-
-
- /**
- * Ask the manifold descriptor to return intermediate points on the object
- * pointed to by the TriaIterator @p iter. This function tries to be
- * backward compatible with respect to the differences between
- * Boundary<dim,spacedim> and Manifold<dim,spacedim>, querying the first
- * whenever the passed @p manifold can be upgraded to a
- * Boundary<dim,spacedim>.
- */
- template <class TriaIterator>
- void get_intermediate_points_on_object(const Manifold<dim, spacedim> &manifold,
- const TriaIterator &iter,
- std::vector<Point<spacedim> > &points) const;
-
/**
* @name Interface with FEValues
* @{
* This function is needed by the constructor of
* <tt>MappingQ<dim,spacedim></tt> for <tt>dim=</tt> 2 and 3.
*
- * For <tt>degree<4</tt> this function sets the @p laplace_on_quad_vector to
+ * For <tt>degree<4</tt> this function sets the @p support_point_weights_on_quad to
* the hardcoded data. For <tt>degree>=4</tt> and MappingQ<2> this vector is
* computed.
*
- * For the definition of the @p laplace_on_quad_vector please refer to
+ * For the definition of the @p support_point_weights_on_quad please refer to
* equation (8) of the `mapping' report.
*/
void
- set_laplace_on_quad_vector(Table<2,double> &loqvs) const;
+ set_support_point_weights_on_quad(Table<2,double> &loqvs) const;
/**
* This function is needed by the constructor of <tt>MappingQ<3></tt>.
*
- * For <tt>degree==2</tt> this function sets the @p laplace_on_hex_vector to
+ * For <tt>degree==2</tt> this function sets the @p support_point_weights_on_hex to
* the hardcoded data. For <tt>degree>2</tt> this vector is computed.
*
- * For the definition of the @p laplace_on_hex_vector please refer to
+ * For the definition of the @p support_point_weights_on_hex please refer to
* equation (8) of the `mapping' report.
*/
- void set_laplace_on_hex_vector(Table<2,double> &lohvs) const;
+ void set_support_point_weights_on_hex(Table<2,double> &lohvs) const;
/**
- * Computes the <tt>laplace_on_quad(hex)_vector</tt>.
+ * Compute the <tt>support_point_weights_on_quad(hex)_vector</tt>.
*
- * Called by the <tt>set_laplace_on_quad(hex)_vector</tt> functions if the
+ * Called by the <tt>set_support_point_weights_on_quad(hex)_vector</tt> functions if the
* data is not yet hardcoded.
*
- * For the definition of the <tt>laplace_on_quad(hex)_vector</tt> please
+ * For the definition of the <tt>support_point_weights_on_quad(hex)_vector</tt> please
* refer to equation (8) of the `mapping' report.
*/
void compute_laplace_vector(Table<2,double> &lvs) const;
/**
- * Takes a <tt>laplace_on_hex(quad)_vector</tt> and applies it to the vector
- * @p a to compute the inner support points as a linear combination of the
- * exterior points.
- *
- * The vector @p a initially contains the locations of the @p n_outer
- * points, the @p n_inner computed inner points are appended.
- *
- * See equation (7) of the `mapping' report.
+ * Compute the support points of the mapping. Interior support
+ * points (ie. support points in quads for 2d, in hexes for 3d) are
+ * computed using the solution of a Laplace equation with the
+ * position of the outer support points as boundary values, in order
+ * to make the transformation as smooth as possible.
*/
- void apply_laplace_vector(const Table<2,double> &lvs,
- std::vector<Point<spacedim> > &a) const;
+ virtual
+ void
+ compute_mapping_support_points(const typename Triangulation<dim,spacedim>::cell_iterator &cell,
+ std::vector<Point<spacedim> > &a) const;
+
/**
- * Computes the support points of the mapping.
- */
- virtual void compute_mapping_support_points(
- const typename Triangulation<dim,spacedim>::cell_iterator &cell,
- std::vector<Point<spacedim> > &a) const;
+ * For <tt>dim=2,3</tt>. Append the support points of all shape functions
+ * located on bounding lines to the vector @p a. Points located on the line
+ * but not on vertices are not included.
+ *
+ * Needed by the @p compute_support_points_laplace function . For
+ * <tt>dim=1</tt> 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.
+ */
+ virtual
+ void
+ add_line_support_points (const typename Triangulation<dim,spacedim>::cell_iterator &cell,
+ std::vector<Point<spacedim> > &a) const;
/**
- * Computes all support points of the mapping shape functions. The inner
- * support points (ie. support points in quads for 2d, in hexes for 3d) are
- * computed using the solution of a Laplace equation with the position of
- * the outer support points as boundary values, in order to make the
- * transformation as smooth as possible.
+ * For <tt>dim=3</tt>. Append the support points of all shape functions
+ * located on bounding faces (quads in 3d) to the vector @p a. Points
+ * located on the quad but not on vertices are not included.
+ *
+ * Needed by the @p compute_support_points_laplace function. For
+ * <tt>dim=1</tt> and <tt>dim=2</tt> 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 compute_support_points_laplace(
- const typename Triangulation<dim,spacedim>::cell_iterator &cell,
- std::vector<Point<spacedim> > &a) const;
+ virtual
+ void
+ add_quad_support_points(const typename Triangulation<dim,spacedim>::cell_iterator &cell,
+ std::vector<Point<spacedim> > &a) const;
/**
- * Needed by the @p laplace_on_quad function (for <tt>dim==2</tt>). Filled
+ * Needed by the @p support_point_weights_on_quad function (for <tt>dim==2</tt>). Filled
* by the constructor.
*
- * Sizes: laplace_on_quad_vector.size()= number of inner unit_support_points
- * laplace_on_quad_vector[i].size()= number of outer unit_support_points,
+ * Sizes: support_point_weights_on_quad.size()= number of inner unit_support_points
+ * support_point_weights_on_quad[i].size()= number of outer unit_support_points,
* i.e. unit_support_points on the boundary of the quad
*
* For the definition of this vector see equation (8) of the `mapping'
* report.
*/
- Table<2,double> laplace_on_quad_vector;
+ Table<2,double> support_point_weights_on_quad;
/**
- * Needed by the @p laplace_on_hex function (for <tt>dim==3</tt>). Filled by
+ * Needed by the @p support_point_weights_on_hex function (for <tt>dim==3</tt>). Filled by
* the constructor.
*
* For the definition of this vector see equation (8) of the `mapping'
* report.
*/
- Table<2,double> laplace_on_hex_vector;
+ Table<2,double> support_point_weights_on_hex;
/**
* Exception.
#ifndef DOXYGEN
template <>
-void MappingQ<1>::set_laplace_on_quad_vector(Table<2,double> &) const;
+void MappingQ<1>::set_support_point_weights_on_quad(Table<2,double> &) const;
template <>
-void MappingQ<3>::set_laplace_on_hex_vector(Table<2,double> &lohvs) const;
+void MappingQ<3>::set_support_point_weights_on_hex(Table<2,double> &lohvs) const;
template <>
void MappingQ<1>::compute_laplace_vector(Table<2,double> &) const;
-template<>
-void MappingQ<3>::add_quad_support_points(const Triangulation<3>::cell_iterator &cell,
- std::vector<Point<3> > &a) const;
-
-// ---- Templated functions ---- //
-template <int dim, int spacedim>
-template <class TriaIterator>
-void
-MappingQ<dim,spacedim>::get_intermediate_points_on_object(const Manifold<dim, spacedim> &manifold,
- const TriaIterator &iter,
- std::vector<Point<spacedim> > &points) const
-{
- const unsigned int structdim = TriaIterator::AccessorType::structure_dimension;
-
- // Try backward compatibility option.
- if (const Boundary<dim,spacedim> *boundary
- = dynamic_cast<const Boundary<dim,spacedim> *>(&manifold))
- // This is actually a boundary. Call old methods.
- {
- switch (structdim)
- {
- case 1:
- {
- const typename Triangulation<dim,spacedim>::line_iterator line = iter;
- boundary->get_intermediate_points_on_line(line, points);
- return;
- }
- case 2:
- {
- const typename Triangulation<dim,spacedim>::quad_iterator quad = iter;
- boundary->get_intermediate_points_on_quad(quad, points);
- return;
- }
- default:
- Assert(false, ExcInternalError());
- return;
- }
- }
- else
- {
- std::vector<Point<spacedim> > sp(GeometryInfo<structdim>::vertices_per_cell);
- for (unsigned int i=0; i<sp.size(); ++i)
- sp[i] = iter->vertex(i);
- get_intermediate_points(manifold, sp, points);
- }
-}
-
-
#endif // DOXYGEN
DEAL_II_NAMESPACE_CLOSE
Assert(n_inner+n_outer==Utilities::fixed_power<dim>(degree+1),
ExcInternalError());
- // build laplace_on_quad_vector
+ // build support_point_weights_on_quad
if (degree>1)
{
if (dim >= 2)
- set_laplace_on_quad_vector(laplace_on_quad_vector);
+ set_support_point_weights_on_quad(support_point_weights_on_quad);
if (dim >= 3)
- set_laplace_on_hex_vector(laplace_on_hex_vector);
+ set_support_point_weights_on_hex(support_point_weights_on_hex);
}
}
Assert(n_inner+n_outer==Utilities::fixed_power<dim>(this->polynomial_degree+1),
ExcInternalError());
- // build laplace_on_quad_vector
+ // build support_point_weights_on_quad
if (this->polynomial_degree>1)
{
if (dim >= 2)
- set_laplace_on_quad_vector(laplace_on_quad_vector);
+ set_support_point_weights_on_quad(support_point_weights_on_quad);
if (dim >= 3)
- set_laplace_on_hex_vector(laplace_on_hex_vector);
+ set_support_point_weights_on_hex(support_point_weights_on_hex);
}
}
template <>
void
-MappingQ<1>::set_laplace_on_quad_vector(Table<2,double> &) const
+MappingQ<1>::set_support_point_weights_on_quad(Table<2,double> &) const
{
Assert(false, ExcInternalError());
}
template<int dim, int spacedim>
void
-MappingQ<dim,spacedim>::set_laplace_on_quad_vector(Table<2,double> &loqvs) const
+MappingQ<dim,spacedim>::set_support_point_weights_on_quad(Table<2,double> &loqvs) const
{
Assert(this->polynomial_degree>1, ExcInternalError());
const unsigned int n_inner_2d=(this->polynomial_degree-1)*(this->polynomial_degree-1);
else if (dim == 3)
{
MappingQ<2,2> mapping_2d(this->polynomial_degree);
- loqvs = mapping_2d.laplace_on_quad_vector;
+ loqvs = mapping_2d.support_point_weights_on_quad;
}
}
template <>
void
-MappingQ<3>::set_laplace_on_hex_vector(Table<2,double> &lohvs) const
+MappingQ<3>::set_support_point_weights_on_hex(Table<2,double> &lohvs) const
{
Assert(this->polynomial_degree>1, ExcInternalError());
template<int dim, int spacedim>
void
-MappingQ<dim,spacedim>::set_laplace_on_hex_vector(Table<2,double> &) const
+MappingQ<dim,spacedim>::set_support_point_weights_on_hex(Table<2,double> &) const
{
Assert(false, ExcInternalError());
}
-template<int dim, int spacedim>
-void
-MappingQ<dim,spacedim>::apply_laplace_vector(const Table<2,double> &lvs,
- std::vector<Point<spacedim> > &a) const
-{
- // check whether the data we need is really available. if you fail here and
- // if lvs==laplace_on_quad_vector in the calling function, then we didn't
- // compute the quad laplace vector. this is mentioned in the constructor of
- // this class, although I don't understand the reason for not aborting there
- // any more [WB]
- Assert(lvs.n_rows()!=0, ExcLaplaceVectorNotSet(this->polynomial_degree));
-
- const unsigned int n_inner_apply=lvs.n_rows();
- Assert(n_inner_apply==n_inner || n_inner_apply==(this->polynomial_degree-1)*(this->polynomial_degree-1),
- ExcInternalError());
- const unsigned int n_outer_apply=lvs.n_cols();
- 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)
- {
- Assert(lvs.n_cols()==n_outer_apply, ExcInternalError());
- Point<spacedim> p;
- for (unsigned int k=0; k<n_outer_apply; ++k)
- p+=lvs[unit_point][k]*a[k];
-
- a.push_back(p);
- }
-}
-
-
-template<int dim, int spacedim>
-void
-MappingQ<dim,spacedim>::compute_mapping_support_points(
- const typename Triangulation<dim,spacedim>::cell_iterator &cell,
- std::vector<Point<spacedim> > &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);
- 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);
-
- for (unsigned int i=0; i<GeometryInfo<dim>::vertices_per_cell; ++i)
- a[i] = cell->vertex(i);
- }
-}
-
-
-template<int dim, int spacedim>
-void
-MappingQ<dim,spacedim>::compute_support_points_laplace(const typename Triangulation<dim,spacedim>::cell_iterator &cell,
- std::vector<Point<spacedim> > &a) const
+namespace
{
- // in any case, we need the vertices first
- a.resize(GeometryInfo<dim>::vertices_per_cell);
- for (unsigned int i=0; i<GeometryInfo<dim>::vertices_per_cell; ++i)
- a[i] = cell->vertex(i);
-
- if (this->polynomial_degree>1)
- switch (dim)
+ /**
+ * Ask the manifold descriptor to return intermediate points on lines or
+ * faces. The function needs to return one or multiple points (depending on
+ * the number of elements in the output vector @p points that lie inside a
+ * line, quad or hex). Whether it is a line, quad or hex doesn't really
+ * matter to this function but it can be inferred from the number of input
+ * points in the @p surrounding_points vector.
+ */
+ template<int dim, int spacedim>
+ void
+ get_intermediate_points (const Manifold<dim, spacedim> &manifold,
+ const QGaussLobatto<1> &line_support_points,
+ const std::vector<Point<spacedim> > &surrounding_points,
+ std::vector<Point<spacedim> > &points)
+ {
+ Assert(surrounding_points.size() >= 2, ExcMessage("At least 2 surrounding points are required"));
+ const unsigned int n=points.size();
+ Assert(n>0, ExcMessage("You can't ask for 0 intermediate points."));
+ std::vector<double> w(surrounding_points.size());
+
+ switch (surrounding_points.size())
{
- case 1:
- add_line_support_points(cell, a);
- break;
case 2:
- // in 2d, add the points on the four bounding lines to the exterior
- // (outer) points
- add_line_support_points (cell, a);
- if (dim != spacedim)
- add_quad_support_points(cell, a);
- else
- apply_laplace_vector (laplace_on_quad_vector,a);
+ {
+ // If two points are passed, these are the two vertices, and
+ // we can only compute degree-1 intermediate points.
+ for (unsigned int i=0; i<n; ++i)
+ {
+ const double x = line_support_points.point(i+1)[0];
+ w[1] = x;
+ w[0] = (1-x);
+ Quadrature<spacedim> quadrature(surrounding_points, w);
+ points[i] = manifold.get_new_point(quadrature);
+ }
break;
+ }
- case 3:
+ case 4:
{
- // in 3d also add the points located on the boundary faces
- add_line_support_points (cell, a);
- add_quad_support_points (cell, a);
- apply_laplace_vector (laplace_on_hex_vector, a);
+ Assert(spacedim >= 2, ExcImpossibleInDim(spacedim));
+ const unsigned m=
+ static_cast<unsigned int>(std::sqrt(static_cast<double>(n)));
+ // is n a square number
+ Assert(m*m==n, ExcInternalError());
+
+ // If four points are passed, these are the two vertices, and
+ // we can only compute (degree-1)*(degree-1) intermediate
+ // points.
+ for (unsigned int i=0; i<m; ++i)
+ {
+ const double y=line_support_points.point(1+i)[0];
+ for (unsigned int j=0; j<m; ++j)
+ {
+ const double x=line_support_points.point(1+j)[0];
+
+ w[0] = (1-x)*(1-y);
+ w[1] = x*(1-y);
+ w[2] = (1-x)*y ;
+ w[3] = x*y ;
+ Quadrature<spacedim> quadrature(surrounding_points, w);
+ points[i*m+j]=manifold.get_new_point(quadrature);
+ }
+ }
break;
}
- default:
+
+ case 8:
Assert(false, ExcNotImplemented());
break;
+ default:
+ Assert(false, ExcInternalError());
+ break;
+ }
+ }
+
+
+
+
+ /**
+ * Ask the manifold descriptor to return intermediate points on the object
+ * pointed to by the TriaIterator @p iter. This function tries to be
+ * backward compatible with respect to the differences between
+ * Boundary<dim,spacedim> and Manifold<dim,spacedim>, querying the first
+ * whenever the passed @p manifold can be upgraded to a
+ * Boundary<dim,spacedim>.
+ */
+ template <int dim, int spacedim, class TriaIterator>
+ void get_intermediate_points_on_object(const Manifold<dim, spacedim> &manifold,
+ const QGaussLobatto<1> &line_support_points,
+ const TriaIterator &iter,
+ std::vector<Point<spacedim> > &points)
+ {
+ const unsigned int structdim = TriaIterator::AccessorType::structure_dimension;
+
+ // Try backward compatibility option.
+ if (const Boundary<dim,spacedim> *boundary
+ = dynamic_cast<const Boundary<dim,spacedim> *>(&manifold))
+ // This is actually a boundary. Call old methods.
+ {
+ switch (structdim)
+ {
+ case 1:
+ {
+ const typename Triangulation<dim,spacedim>::line_iterator line = iter;
+ boundary->get_intermediate_points_on_line(line, points);
+ return;
+ }
+ case 2:
+ {
+ const typename Triangulation<dim,spacedim>::quad_iterator quad = iter;
+ boundary->get_intermediate_points_on_quad(quad, points);
+ return;
+ }
+ default:
+ Assert(false, ExcInternalError());
+ return;
+ }
}
+ else
+ {
+ std::vector<Point<spacedim> > sp(GeometryInfo<structdim>::vertices_per_cell);
+ for (unsigned int i=0; i<sp.size(); ++i)
+ sp[i] = iter->vertex(i);
+ get_intermediate_points(manifold, line_support_points, sp, points);
+ }
+ }
+
+
+ /**
+ * Take a <tt>support_point_weights_on_hex(quad)</tt> and apply it to the vector
+ * @p a to compute the inner support points as a linear combination of the
+ * exterior points.
+ *
+ * The vector @p a initially contains the locations of the @p n_outer
+ * points, the @p n_inner computed inner points are appended.
+ *
+ * See equation (7) of the `mapping' report.
+ */
+ template <int spacedim>
+ void add_weighted_interior_points(const Table<2,double> &lvs,
+ std::vector<Point<spacedim> > &a)
+ {
+ const unsigned int n_inner_apply=lvs.n_rows();
+ const unsigned int n_outer_apply=lvs.n_cols();
+ 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)
+ {
+ Assert(lvs.n_cols()==n_outer_apply, ExcInternalError());
+ Point<spacedim> p;
+ for (unsigned int k=0; k<n_outer_apply; ++k)
+ p+=lvs[unit_point][k]*a[k];
+
+ a.push_back(p);
+ }
+ }
}
-template<int dim, int spacedim>
+template <int dim, int spacedim>
void
-MappingQ<dim,spacedim>::add_line_support_points (const typename Triangulation<dim,spacedim>::cell_iterator &cell,
- std::vector<Point<spacedim> > &a) const
+MappingQ<dim,spacedim>::
+add_line_support_points (const typename Triangulation<dim,spacedim>::cell_iterator &cell,
+ std::vector<Point<spacedim> > &a) const
{
// if we only need the midpoint, then ask for it.
if (this->polynomial_degree==2)
const Manifold<dim,spacedim> &manifold =
( ( line->manifold_id() == numbers::invalid_manifold_id ) &&
- ( dim < spacedim ) ? cell->get_manifold() :
+ ( dim < spacedim )
+ ?
+ cell->get_manifold()
+ :
line->get_manifold() );
a.push_back(manifold.get_new_point_on_line(line));
- };
+ }
}
else
// otherwise call the more complicated functions and ask for inner points
cell->get_manifold() :
line->get_manifold() );
- get_intermediate_points_on_object (manifold, line, line_points);
+ get_intermediate_points_on_object (manifold, line_support_points, line, line_points);
if (dim==3)
{
// in 2D, lines always have the correct orientation. simply append
// all points
a.insert (a.end(), line_points.begin(), line_points.end());
-
}
}
}
-template<>
+template <>
void
-MappingQ<3>::
+MappingQ<3,3>::
add_quad_support_points(const Triangulation<3>::cell_iterator &cell,
std::vector<Point<3> > &a) const
{
static const StraightBoundary<3> straight_boundary;
// used if face quad at boundary or entirely in the interior of the domain
- std::vector<Point<3> > quad_points ((this->polynomial_degree-1)*(this->polynomial_degree-1));
+ std::vector<Point<3> > quad_points ((polynomial_degree-1)*(polynomial_degree-1));
// used if only one line of face quad is at boundary
- std::vector<Point<3> > b(4*this->polynomial_degree);
+ std::vector<Point<3> > b(4*polynomial_degree);
// Used by the new Manifold interface. This vector collects the
// vertices used to compute the intermediate points.
// points on it
if (face->at_boundary())
{
- get_intermediate_points_on_object(face->get_manifold(), face, quad_points);
+ get_intermediate_points_on_object(face->get_manifold(), line_support_points, face, quad_points);
// in 3D, the orientation, flip and rotation of the face might not
// match what we expect here, namely the standard orientation. thus
// reorder points accordingly. since a Mapping uses the same shape
- // function as an FEQ, we can ask a FEQ to do the reordering for us.
+ // function as an FE_Q, we can ask a FE_Q to do the reordering for us.
for (unsigned int i=0; i<quad_points.size(); ++i)
a.push_back(quad_points[feq.adjust_quad_dof_index_for_face_orientation(i,
face_orientation,
// boundary, then collect points separately
if (lines_at_boundary>0)
{
- // call of function apply_laplace_vector increases size of b
+ // call of function add_weighted_interior_points increases size of b
// about 1. There resize b for the case the mentioned function
// was already called.
- b.resize(4*this->polynomial_degree);
+ b.resize(4*polynomial_degree);
// b is of size 4*degree, make sure that this is the right size
- Assert(b.size()==vertices_per_face+lines_per_face*(this->polynomial_degree-1),
+ Assert(b.size()==vertices_per_face+lines_per_face*(polynomial_degree-1),
ExcDimensionMismatch(b.size(),
- vertices_per_face+lines_per_face*(this->polynomial_degree-1)));
+ vertices_per_face+lines_per_face*(polynomial_degree-1)));
// sort the points into b. We used access from the cell (not
// from the face) to fill b, so we can assume a standard face
b[i]=a[GeometryInfo<3>::face_to_cell_vertices(face_no, i)];
for (unsigned int i=0; i<lines_per_face; ++i)
- for (unsigned int j=0; j<this->polynomial_degree-1; ++j)
- b[vertices_per_face+i*(this->polynomial_degree-1)+j]=
+ for (unsigned int j=0; j<polynomial_degree-1; ++j)
+ b[vertices_per_face+i*(polynomial_degree-1)+j]=
a[vertices_per_cell + GeometryInfo<3>::face_to_cell_lines(
- face_no, i)*(this->polynomial_degree-1)+j];
+ face_no, i)*(polynomial_degree-1)+j];
// Now b includes the support points on the quad and we can
// apply the laplace vector
- apply_laplace_vector(laplace_on_quad_vector, b);
- Assert(b.size()==4*this->polynomial_degree+(this->polynomial_degree-1)*(this->polynomial_degree-1),
- ExcDimensionMismatch(b.size(), 4*this->polynomial_degree+(this->polynomial_degree-1)*(this->polynomial_degree-1)));
+ add_weighted_interior_points (support_point_weights_on_quad, b);
+ AssertDimension (b.size(),
+ 4*this->polynomial_degree +
+ (this->polynomial_degree-1)*(this->polynomial_degree-1));
- for (unsigned int i=0; i<(this->polynomial_degree-1)*(this->polynomial_degree-1); ++i)
- a.push_back(b[4*this->polynomial_degree+i]);
+ for (unsigned int i=0; i<(polynomial_degree-1)*(polynomial_degree-1); ++i)
+ a.push_back(b[4*polynomial_degree+i]);
}
else
{
vertices.resize(4);
for (unsigned int i=0; i<4; ++i)
vertices[i] = face->vertex(i);
- get_intermediate_points (face->get_manifold(), vertices, quad_points);
+ get_intermediate_points (face->get_manifold(), line_support_points, vertices, quad_points);
// in 3D, the orientation, flip and rotation of the face might
// not match what we expect here, namely the standard
// orientation. thus reorder points accordingly. since a Mapping
-template<>
+template <>
void
MappingQ<2,3>::
add_quad_support_points(const Triangulation<2,3>::cell_iterator &cell,
std::vector<Point<3> > &a) const
{
- std::vector<Point<3> > quad_points ((this->polynomial_degree-1)*(this->polynomial_degree-1));
- get_intermediate_points_on_object (cell->get_manifold(), cell, quad_points);
+ std::vector<Point<3> > quad_points ((polynomial_degree-1)*(polynomial_degree-1));
+ get_intermediate_points_on_object (cell->get_manifold(), line_support_points,
+ cell, quad_points);
for (unsigned int i=0; i<quad_points.size(); ++i)
a.push_back(quad_points[i]);
}
-template<int dim, int spacedim>
+template <int dim, int spacedim>
void
MappingQ<dim,spacedim>::
add_quad_support_points(const typename Triangulation<dim,spacedim>::cell_iterator &,
std::vector<Point<spacedim> > &) const
{
- Assert (dim > 2, ExcImpossibleInDim(dim));
+ Assert (false, ExcInternalError());
+}
+
+
+
+
+
+template<int dim, int spacedim>
+void
+MappingQ<dim,spacedim>::
+compute_mapping_support_points(const typename Triangulation<dim,spacedim>::cell_iterator &cell,
+ std::vector<Point<spacedim> > &a) const
+{
+ // get the vertices first
+ a.resize(GeometryInfo<dim>::vertices_per_cell);
+ for (unsigned int i=0; i<GeometryInfo<dim>::vertices_per_cell; ++i)
+ a[i] = cell->vertex(i);
+
+ if (this->polynomial_degree>1)
+ switch (dim)
+ {
+ case 1:
+ add_line_support_points(cell, a);
+ break;
+ case 2:
+ // in 2d, add the points on the four bounding lines to the exterior
+ // (outer) points
+ add_line_support_points(cell, a);
+
+ // then get the support points on the quad if we are on a
+ // manifold, otherwise compute them from the points around it
+ if (dim != spacedim)
+ add_quad_support_points(cell, a);
+ else
+ add_weighted_interior_points (support_point_weights_on_quad, a);
+ break;
+
+ case 3:
+ {
+ // in 3d also add the points located on the boundary faces
+ add_line_support_points (cell, a);
+ add_quad_support_points (cell, a);
+
+ // then compute the interior points
+ add_weighted_interior_points (support_point_weights_on_hex, a);
+ break;
+ }
+
+ default:
+ Assert(false, ExcNotImplemented());
+ break;
+ }
}
}
-template<int dim, int spacedim>
-void
-MappingQ<dim,spacedim>::get_intermediate_points (const Manifold<dim, spacedim> &manifold,
- const std::vector<Point<spacedim> > &surrounding_points,
- std::vector<Point<spacedim> > &points) const
-{
- Assert(surrounding_points.size() >= 2, ExcMessage("At least 2 surrounding points are required"));
- const unsigned int n=points.size();
- Assert(n>0, ExcMessage("You can't ask for 0 intermediate points."));
- std::vector<double> w(surrounding_points.size());
-
- switch (surrounding_points.size())
- {
- case 2:
- {
- // If two points are passed, these are the two vertices, and
- // we can only compute degree-1 intermediate points.
- AssertDimension(n, this->polynomial_degree-1);
- for (unsigned int i=0; i<n; ++i)
- {
- const double x = line_support_points.point(i+1)[0];
- w[1] = x;
- w[0] = (1-x);
- Quadrature<spacedim> quadrature(surrounding_points, w);
- points[i] = manifold.get_new_point(quadrature);
- }
- break;
- }
-
- case 4:
- {
- Assert(spacedim >= 2, ExcImpossibleInDim(spacedim));
- const unsigned m=
- static_cast<unsigned int>(std::sqrt(static_cast<double>(n)));
- // is n a square number
- Assert(m*m==n, ExcInternalError());
-
- // If four points are passed, these are the two vertices, and
- // we can only compute (degree-1)*(degree-1) intermediate
- // points.
- AssertDimension(m, this->polynomial_degree-1);
-
- for (unsigned int i=0; i<m; ++i)
- {
- const double y=line_support_points.point(1+i)[0];
- for (unsigned int j=0; j<m; ++j)
- {
- const double x=line_support_points.point(1+j)[0];
-
- w[0] = (1-x)*(1-y);
- w[1] = x*(1-y);
- w[2] = (1-x)*y ;
- w[3] = x*y ;
- Quadrature<spacedim> quadrature(surrounding_points, w);
- points[i*m+j]=manifold.get_new_point(quadrature);
- }
- }
- break;
- }
-
- case 8:
- Assert(false, ExcNotImplemented());
- break;
- default:
- Assert(false, ExcInternalError());
- break;
- }
-}
-
-
// explicit instantiations
#include "mapping_q.inst"
projects / dealii.git / commitdiff