*
* To use this class, make sure that the
* Boundary::@p get_normals_at_vertices function is implemented
- * for the users boundary object.
+ * for the user's boundary object.
*
* For more information about the <tt>spacedim</tt> template parameter
* check the documentation of FiniteElement or the one of
{
public:
/**
- * Constructor. Pass the fixed
- * degree @p 3 down to the base
- * class, as a cubic mapping
- * suffices to generate a
- * continuous mapping of the
- * boundary.
+ * Constructor. Pass the fixed degree @p 3 down to the base class, as a
+ * cubic mapping suffices to generate a continuous mapping of the boundary.
*/
MappingC1 ();
/**
- * Return a pointer to a copy of the
- * present object. The caller of this
- * copy then assumes ownership of it.
+ * Return a pointer to a copy of the present object. The caller of this copy
+ * then assumes ownership of it.
*/
virtual
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 on vertices are not
- * included.
+ * 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 on vertices are not included.
*
- * Needed by the
- * <tt>compute_support_points_simple(laplace)</tt>
- * functions. For <tt>dim=1</tt> this
- * function is empty.
+ * Needed by the <tt>compute_support_points_simple(laplace)</tt>
+ * functions. For <tt>dim=1</tt> this function is empty.
*
- * This function chooses the
- * respective points not such
- * that they are interpolating
- * the boundary (as does the base
- * class), but rather such that
- * the resulting cubic mapping is
- * a continuous one.
+ * This function chooses the respective points not such that they are
+ * interpolating the boundary (as does the base class), but rather such that
+ * the resulting cubic mapping is a continuous one.
*/
virtual void
add_line_support_points (const typename Triangulation<dim>::cell_iterator &cell,
std::vector<Point<dim> > &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 line but on vertices
- * are not included.
+ * 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 line but on vertices are not included.
*
- * Needed by the
- * @p compute_support_points_laplace
- * function. For <tt>dim=1</tt> and 2
- * this function is empty.
+ * Needed by the @p compute_support_points_laplace function. For
+ * <tt>dim=1</tt> and 2 this function is empty.
*
- * This function chooses the
- * respective points not such
- * that they are interpolating
- * the boundary (as does the base
- * class), but rather such that
- * the resulting cubic mapping is
- * a continuous one.
+ * This function chooses the respective points not such that they are
+ * interpolating the boundary (as does the base class), but rather such that
+ * the resulting cubic mapping is a continuous one.
*/
virtual void
add_quad_support_points(const typename Triangulation<dim>::cell_iterator &cell,
/*@{*/
/**
- * Mapping class that uses Qp-mappings on boundary cells. The mapping
- * shape functions make use of tensor product polynomials with
- * equidistant (on the unit cell) support points.
+ * Mapping class that uses Qp-mappings on boundary cells. The mapping shape
+ * functions make use of tensor product polynomials with unit cell support
+ * points equal to the points of the Gauss-Lobatto quadrature formula. These
+ * points give a well-conditioned interpolation also for very high orders and
+ * are therefore preferred over equidistant support points.
*
* For more details about Qp-mappings, see the `mapping' report at
* <tt>deal.II/doc/reports/mapping_q/index.html</tt> in the `Reports'
* check the documentation of FiniteElement or the one of
* Triangulation.
*
+ * @note Since the boundary description is closely tied to the unit cell
+ * support points, new boundary descriptions need to explicitly use the
+ * Gauss-Lobatto points.
+ *
* @author Ralf Hartmann, 2000, 2001, 2005; Guido Kanschat 2000, 2001
*/
template <int dim, int spacedim=dim>
{
public:
/**
- * Constructor. @p p gives the
- * degree of mapping polynomials
- * on boundary cells.
+ * Constructor. @p p gives the degree of mapping polynomials on boundary
+ * cells.
*
- * The second argument determines
- * whether the higher order
- * mapping should also be used on
- * interior cells. If its value
- * is <code>false</code> (the
- * default), the a lower-order
- * mapping is used in the
- * interior. This is sufficient
- * for most cases where higher
- * order mappings are only used
- * to better approximate the
- * boundary. In that case, cells
- * bounded by straight lines are
- * acceptable in the
- * interior. However, there are
- * cases where one would also
- * like to use a higher order
- * mapping in the interior. The
- * MappingQEulerian class is one
- * such case.
+ * The second argument determines whether the higher order mapping should
+ * also be used on interior cells. If its value is <code>false</code> (the
+ * default), the a lower-order mapping is used in the interior. This is
+ * sufficient for most cases where higher order mappings are only used to
+ * better approximate the boundary. In that case, cells bounded by straight
+ * lines are acceptable in the interior. However, there are cases where one
+ * would also like to use a higher order mapping in the interior. The
+ * MappingQEulerian class is one such case.
*/
MappingQ (const unsigned int p,
const bool use_mapping_q_on_all_cells = false);
/**
- * Copy constructor. Performs a
- * deep copy, i.e. duplicates
- * what #tensor_pols points to
- * instead of simply copying the
- * #tensor_pols pointer as done
- * by a default copy constructor.
+ * Copy constructor. Performs a deep copy, i.e. duplicates what #tensor_pols
+ * points to instead of simply copying the #tensor_pols pointer as done by a
+ * default copy constructor.
*/
MappingQ (const MappingQ<dim,spacedim> &mapping);
virtual ~MappingQ ();
/**
- * Transforms the point @p p on
- * the unit cell to the point
- * @p p_real on the real cell
- * @p cell and returns @p p_real.
+ * Transforms the point @p p on the unit cell to the point @p p_real on the
+ * real cell @p cell and returns @p p_real.
*/
virtual Point<spacedim>
transform_unit_to_real_cell (
const Point<dim> &p) const;
/**
- * Transforms the point @p p on
- * the real cell to the point
- * @p p_unit on the unit cell
- * @p cell and returns @p p_unit.
+ * Transforms the point @p p on the real cell to the point @p p_unit on the
+ * unit cell @p cell and returns @p p_unit.
*
- * Uses Newton iteration and the
- * @p transform_unit_to_real_cell
- * function.
+ * Uses Newton iteration and the @p transform_unit_to_real_cell function.
*
- * In the codimension one case,
- * this function returns the
- * normal projection of the real
- * point @p p on the curve or
- * surface identified by the @p
- * cell.
+ * In the codimension one case, this function returns the normal projection
+ * of the real point @p p on the curve or surface identified by the @p cell.
*
- * @note Polynomial mappings from
- * the reference (unit) cell coordinates
- * to the coordinate system of a real
- * cell are not always invertible if
- * the point for which the inverse
- * mapping is to be computed lies
- * outside the cell's boundaries.
- * In such cases, the current function
- * may fail to compute a point on
- * the reference cell whose image
- * under the mapping equals the given
- * point @p p. If this is the case
- * then this function throws an
- * exception of type
- * Mapping::ExcTransformationFailed .
- * Whether the given point @p p lies
- * outside the cell can therefore be
- * determined by checking whether the
- * return reference coordinates lie
- * inside of outside the reference
- * cell (e.g., using
- * GeometryInfo::is_inside_unit_cell)
- * or whether the exception mentioned
- * above has been thrown.
+ * @note Polynomial mappings from the reference (unit) cell coordinates to
+ * the coordinate system of a real cell are not always invertible if the
+ * point for which the inverse mapping is to be computed lies outside the
+ * cell's boundaries. In such cases, the current function may fail to
+ * compute a point on the reference cell whose image under the mapping
+ * equals the given point @p p. If this is the case then this function
+ * throws an exception of type Mapping::ExcTransformationFailed . Whether
+ * the given point @p p lies outside the cell can therefore be determined by
+ * checking whether the return reference coordinates lie inside of outside
+ * the reference cell (e.g., using GeometryInfo::is_inside_unit_cell) or
+ * whether the exception mentioned above has been thrown.
*/
virtual Point<dim>
transform_real_to_unit_cell (
const MappingType type) const;
/**
- * Return the degree of the
- * mapping, i.e. the value which
- * was passed to the constructor.
+ * Return the degree of the mapping, i.e. the value which was passed to the
+ * constructor.
*/
unsigned int get_degree () const;
/**
- * Return a pointer to a copy of the
- * present object. The caller of this
- * copy then assumes ownership of it.
+ * Return a pointer to a copy of the present object. The caller of this copy
+ * then assumes ownership of it.
*/
virtual
Mapping<dim,spacedim> *clone () const;
/**
- * Storage for internal data of
- * Q_degree transformation.
+ * Storage for internal data of Q_degree transformation.
*/
class InternalData : public MappingQ1<dim,spacedim>::InternalData
{
/**
- * Return an estimate (in
- * bytes) or the memory
- * consumption of this
- * object.
+ * Return an estimate (in bytes) or the memory consumption of this object.
*/
virtual std::size_t memory_consumption () const;
/**
- * Unit normal vectors. Used
- * for the alternative
- * computation of the normal
- * vectors. See doc of the
- * @p alternative_normals_computation
- * flag.
+ * Unit normal vectors. Used for the alternative computation of the normal
+ * vectors. See doc of the @p alternative_normals_computation flag.
*
- * Filled (hardcoded) once in
- * @p get_face_data.
+ * Filled (hardcoded) once in @p get_face_data.
*/
std::vector<std::vector<Point<dim> > > unit_normals;
/**
- * Flag that is set by the
- * <tt>fill_fe_[[sub]face]_values</tt>
- * function.
+ * Flag that is set by the <tt>fill_fe_[[sub]face]_values</tt> function.
*
- * If this flag is @p true
- * we are on an interior cell
- * and the
- * @p mapping_q1_data is
- * used.
+ * If this flag is @p true we are on an interior cell and the @p
+ * mapping_q1_data is used.
*/
bool use_mapping_q1_on_current_cell;
/**
- * On interior cells
- * @p MappingQ1 is used.
+ * On interior cells @p MappingQ1 is used.
*/
typename MappingQ1<dim,spacedim>::InternalData mapping_q1_data;
};
protected:
/**
- * Implementation of the interface in
- * Mapping.
+ * Implementation of the interface in Mapping.
*/
virtual void
fill_fe_values (const typename Triangulation<dim,spacedim>::cell_iterator &cell,
CellSimilarity::Similarity &cell_similarity) const ;
/**
- * Implementation of the interface in
- * Mapping.
+ * Implementation of the interface in Mapping.
*/
virtual void
fill_fe_face_values (const typename Triangulation<dim,spacedim>::cell_iterator &cell,
typename std::vector<Point<spacedim> > &normal_vectors) const ;
/**
- * Implementation of the interface in
- * Mapping.
+ * Implementation of the interface in Mapping.
*/
virtual void
fill_fe_subface_values (const typename Triangulation<dim,spacedim>::cell_iterator &cell,
typename std::vector<Point<spacedim> > &normal_vectors) 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.
+ * 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.
+ * 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.
+ * 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.
+ * 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.
+ * 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.
+ * 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,
const Quadrature<dim-1>& quadrature) const;
/**
- * Compute shape values and/or
- * derivatives.
+ * Compute shape values and/or derivatives.
*/
virtual void
compute_shapes_virtual (const std::vector<Point<dim> > &unit_points,
typename MappingQ1<dim,spacedim>::InternalData &data) const;
/**
- * This function is needed by the
- * constructor of <tt>MappingQ<dim,spacedim></tt>
- * for <tt>dim=</tt> 2 and 3.
+ * 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
- * the hardcoded data. For
- * <tt>degree>=4</tt> and MappingQ<2>
- * this vector is computed.
+ * For <tt>degree<4</tt> this function sets the @p laplace_on_quad_vector 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 equation (8)
- * of the `mapping' report.
+ * For the definition of the @p laplace_on_quad_vector please refer to
+ * equation (8) of the `mapping' report.
*/
void
set_laplace_on_quad_vector(Table<2,double> &loqvs) const;
/**
- * This function is needed by the
- * constructor of <tt>MappingQ<3></tt>.
+ * 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
- * the hardcoded data. For
- * <tt>degree>2</tt> this vector is
- * computed.
+ * For <tt>degree==2</tt> this function sets the @p laplace_on_hex_vector 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 equation (8)
- * of the `mapping' report.
+ * For the definition of the @p laplace_on_hex_vector please refer to
+ * equation (8) of the `mapping' report.
*/
void set_laplace_on_hex_vector(Table<2,double> &lohvs) const;
/**
- * Computes the
- * <tt>laplace_on_quad(hex)_vector</tt>.
+ * Computes the <tt>laplace_on_quad(hex)_vector</tt>.
*
- * Called by the
- * <tt>set_laplace_on_quad(hex)_vector</tt>
- * functions if the data is not
- * yet hardcoded.
+ * Called by the <tt>set_laplace_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 refer to equation (8)
- * of the `mapping' report.
+ * For the definition of the <tt>laplace_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.
+ * 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.
+ * 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.
+ * See equation (7) of the `mapping' report.
*/
void apply_laplace_vector(const Table<2,double> &lvs,
std::vector<Point<spacedim> > &a) const;
/**
- * Computes the support points of
- * the mapping.
+ * 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;
/**
- * 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.
+ * 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.
*/
void compute_support_points_laplace(
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 by the
- * constructor.
+ * Needed by the @p laplace_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, i.e.
- * unit_support_points on the
- * boundary of the quad
+ * laplace_on_quad_vector.size()= number of inner unit_support_points
+ * laplace_on_quad_vector[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.
+ * For the definition of this vector see equation (8) of the `mapping'
+ * report.
*/
Table<2,double> laplace_on_quad_vector;
/**
- * Needed by the
- * @p laplace_on_hex function
- * (for <tt>dim==3</tt>). Filled by the
- * constructor.
+ * Needed by the @p laplace_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.
+ * For the definition of this vector see equation (8) of the `mapping'
+ * report.
*/
Table<2,double> laplace_on_hex_vector;
<< "laplace_vector not set for degree=" << arg1 << ".");
/**
- * Degree @p p of the
- * polynomials used as shape
- * functions for the Qp mapping
+ * Degree @p p of the polynomials used as shape functions for the Qp mapping
* of cells at the boundary.
*/
const unsigned int degree;
/**
- * Number of inner mapping shape
- * functions.
+ * Number of inner mapping shape functions.
*/
const unsigned int n_inner;
/**
- * Number of mapping shape
- * functions on the boundary.
+ * Number of mapping shape functions on the boundary.
*/
const unsigned int n_outer;
/**
- * Pointer to the
- * @p dim-dimensional tensor
- * product polynomials used as
- * shape functions for the Qp
- * mapping of cells at the
- * boundary.
+ * Pointer to the @p dim-dimensional tensor product polynomials used as
+ * shape functions for the Qp mapping of cells at the boundary.
*/
const TensorProductPolynomials<dim> *tensor_pols;
/**
- * Number of the Qp tensor
- * product shape functions.
+ * Number of the Qp tensor product shape functions.
*/
const unsigned int n_shape_functions;
/**
- * Mapping from lexicographic to
- * to the Qp shape function
- * numbering. Its size is
- * @p dofs_per_cell.
+ * Mapping from lexicographic to to the Qp shape function numbering. Its
+ * size is @p dofs_per_cell.
*/
const std::vector<unsigned int> renumber;
/**
- * If this flag is set @p true
- * then @p MappingQ is used on
- * all cells, not only on
- * boundary cells.
+ * If this flag is set @p true then @p MappingQ is used on all cells, not
+ * only on boundary cells.
*/
const bool use_mapping_q_on_all_cells;
/**
- * An FE_Q object which is only needed in
- * 3D, since it knows how to reorder shape
- * functions/DoFs on non-standard
- * faces. This is used to reorder support
- * points in the same way. We could make
- * this a pointer to prevent construction
- * in 1D and 2D, but since memory and time
- * requirements are not particularly high
- * this seems unnecessary at the moment.
+ * An FE_Q object which is only needed in 3D, since it knows how to reorder
+ * shape functions/DoFs on non-standard faces. This is used to reorder
+ * support points in the same way. We could make this a pointer to prevent
+ * construction in 1D and 2D, but since memory and time requirements are not
+ * particularly high this seems unnecessary at the moment.
*/
const FE_Q<dim> feq;
+
+ /**
+ * Declare other MappingQ classes friends.
+ */
+ template <int,int> friend class MappingQ;
};
/*@}*/
{
public:
/**
- * Constructor. The first argument is
- * the polynomical degree of the desired
- * Qp mapping. It then takes a
- * <tt>Vector<double> &</tt> to specify the
- * transformation of the domain
- * from the reference to
- * the current configuration.
- * The organization of the
- * elements in the @p Vector
- * must follow the concept how
- * deal.II stores solutions that
- * are associated to a
- * triangulation. This is
- * automatically the case if the
- * @p Vector represents the
- * solution of the previous step
- * of a nonlinear problem.
- * Alternatively, the @p Vector
- * can be initialized by
+ * Constructor. The first argument is the polynomical degree of the desired
+ * Qp mapping. It then takes a <tt>Vector<double> &</tt> to specify the
+ * transformation of the domain from the reference to the current
+ * configuration. The organization of the elements in the @p Vector must
+ * follow the concept how deal.II stores solutions that are associated to a
+ * triangulation. This is automatically the case if the @p Vector
+ * represents the solution of the previous step of a nonlinear problem.
+ * Alternatively, the @p Vector can be initialized by
* <tt>DoFAccessor::set_dof_values()</tt>.
*/
const DoFHandler<dim,spacedim> &euler_dof_handler);
/**
- * Return a pointer to a copy of the
- * present object. The caller of this
- * copy then assumes ownership of it.
+ * Return a pointer to a copy of the present object. The caller of this copy
+ * then assumes ownership of it.
*/
virtual
Mapping<dim,spacedim> *clone () const;
/**
- * Always returns @p false because
- * MappingQ1Eulerian does not in general
- * preserve vertex locations (unless the
- * translation vector happens to provide
- * for zero displacements at vertex
- * locations).
+ * Always returns @p false because MappingQ1Eulerian does not in general
+ * preserve vertex locations (unless the translation vector happens to
+ * provide for zero displacements at vertex locations).
*/
bool preserves_vertex_locations () const;
protected:
/**
- * Implementation of the interface in
- * MappingQ. Overrides the function in
- * the base class, since we cannot use
- * any cell similarity for this class.
+ * Implementation of the interface in MappingQ. Overrides the function in
+ * the base class, since we cannot use any cell similarity for this class.
*/
virtual void
fill_fe_values (const typename Triangulation<dim,spacedim>::cell_iterator &cell,
CellSimilarity::Similarity &cell_similarity) const;
/**
- * Reference to the vector of
- * shifts.
+ * Reference to the vector of shifts.
*/
SmartPointer<const VECTOR, MappingQEulerian<dim,VECTOR,spacedim> > euler_vector;
/**
- * Pointer to the DoFHandler to
- * which the mapping vector is
- * associated.
+ * Pointer to the DoFHandler to which the mapping vector is associated.
*/
-
SmartPointer<const DoFHandler<dim,spacedim>,MappingQEulerian<dim,VECTOR,spacedim> > euler_dof_handler;
private:
/**
- * Special quadrature rule used
- * to define the support points
- * in the reference configuration.
+ * Special quadrature rule used to define the support points in the
+ * reference configuration.
*/
class SupportQuadrature : public Quadrature<dim>
{
public:
/**
- * Constructor, with an argument
- * defining the desired polynomial
- * degree.
+ * Constructor, with an argument defining the desired polynomial degree.
*/
SupportQuadrature (const unsigned int map_degree);
};
/**
- * A member variable holding the
- * quadrature points in the right
- * order.
+ * A member variable holding the quadrature points in the right order.
*/
const SupportQuadrature support_quadrature;
/**
- * FEValues object used to query the
- * the given finite element field
- * at the support points in the
- * reference configuration.
+ * FEValues object used to query the the given finite element field at the
+ * support points in the reference configuration.
*
- * The variable is marked as
- * mutable since we have to call
- * FEValues::reinit from
- * compute_mapping_support_points,
- * a function that is 'const'.
+ * The variable is marked as mutable since we have to call FEValues::reinit
+ * from compute_mapping_support_points, a function that is 'const'.
*/
mutable FEValues<dim,spacedim> fe_values;
/**
- * A variable to guard access to
- * the fe_values variable.
+ * A variable to guard access to the fe_values variable.
*/
mutable Threads::Mutex fe_values_mutex;
/**
- * Compute the positions of the
- * support points in the current
- * configuration
+ * Compute the positions of the support points in the current configuration
*/
virtual void compute_mapping_support_points(
const typename Triangulation<dim,spacedim>::cell_iterator &cell,
virtual
void
create_notification (const Triangulation<dim, spacedim> &tria);
- } DEAL_II_DEPRECATED;
+ };
/**
* A structure that is used as an
#include <deal.II/base/config.h>
#include <deal.II/base/subscriptor.h>
+#include <deal.II/base/quadrature_lib.h>
+#include <deal.II/base/thread_management.h>
#include <deal.II/base/point.h>
#include <deal.II/grid/tria.h>
public:
/**
- * Type keeping information about
- * the normals at the vertices of
- * a face of a cell. Thus, there
- * are
- * <tt>GeometryInfo<dim>::vertices_per_face</tt>
- * normal vectors, that define
- * the tangent spaces of the
- * boundary at the vertices. Note
- * that the vectors stored in
- * this object are not required
- * to be normalized, nor to
- * actually point outward, as one
- * often will only want to check
- * for orthogonality to define
- * the tangent plane; if a
- * function requires the normals
- * to be normalized, then it must
- * do so itself.
+ * Type keeping information about the normals at the vertices of a face of a
+ * cell. Thus, there are <tt>GeometryInfo<dim>::vertices_per_face</tt>
+ * normal vectors, that define the tangent spaces of the boundary at the
+ * vertices. Note that the vectors stored in this object are not required to
+ * be normalized, nor to actually point outward, as one often will only want
+ * to check for orthogonality to define the tangent plane; if a function
+ * requires the normals to be normalized, then it must do so itself.
*
- * For obvious reasons, this
- * type is not useful in 1d.
+ * For obvious reasons, this type is not useful in 1d.
*/
typedef Tensor<1,spacedim> FaceVertexNormals[GeometryInfo<dim>::vertices_per_face];
/**
- * Destructor. Does nothing here, but
- * needs to be declared to make it
+ * Destructor. Does nothing here, but needs to be declared to make it
* virtual.
*/
virtual ~Boundary ();
/**
- * Return the point which shall
- * become the new middle vertex
- * of the two children of a
- * regular line. In 2D, this line
- * is a line at the boundary,
- * while in 3d, it is bounding a
- * face at the boundary (the
- * lines therefore is also on the
- * boundary).
+ * Return the point which shall become the new middle vertex of the two
+ * children of a regular line. In 2D, this line is a line at the boundary,
+ * while in 3d, it is bounding a face at the boundary (the lines therefore
+ * is also on the boundary).
*/
virtual
Point<spacedim>
get_new_point_on_line (const typename Triangulation<dim,spacedim>::line_iterator &line) const = 0;
/**
- * Return the point which shall
- * become the common point of the
- * four children of a quad at the
- * boundary in three or more
- * spatial dimensions. This
- * function therefore is only
- * useful in at least three
- * dimensions and should not be
- * called for lower dimensions.
+ * Return the point which shall become the common point of the four children
+ * of a quad at the boundary in three or more spatial dimensions. This
+ * function therefore is only useful in at least three dimensions and should
+ * not be called for lower dimensions.
*
- * This function is called after
- * the four lines bounding the
- * given @p quad are refined, so
- * you may want to use the
- * information provided by
- * <tt>quad->line(i)->child(j)</tt>,
- * <tt>i=0...3</tt>, <tt>j=0,1</tt>.
+ * This function is called after the four lines bounding the given @p quad
+ * are refined, so you may want to use the information provided by
+ * <tt>quad->line(i)->child(j)</tt>, <tt>i=0...3</tt>, <tt>j=0,1</tt>.
*
- * Because in 2D, this function
- * is not needed, it is not made
- * pure virtual, to avoid the
- * need to overload it. The
- * default implementation throws
- * an error in any case, however.
+ * Because in 2D, this function is not needed, it is not made pure virtual,
+ * to avoid the need to overload it. The default implementation throws an
+ * error in any case, however.
*/
virtual
Point<spacedim>
get_new_point_on_quad (const typename Triangulation<dim,spacedim>::quad_iterator &quad) const;
/**
- * Depending on <tt>dim=2</tt> or
- * <tt>dim=3</tt> this function
- * calls the
- * get_new_point_on_line or the
- * get_new_point_on_quad
- * function. It throws an
- * exception for
- * <tt>dim=1</tt>. This wrapper
- * allows dimension independent
+ * Depending on <tt>dim=2</tt> or <tt>dim=3</tt> this function calls the
+ * get_new_point_on_line or the get_new_point_on_quad function. It throws an
+ * exception for <tt>dim=1</tt>. This wrapper allows dimension independent
* programming.
*/
Point<spacedim>
get_new_point_on_face (const typename Triangulation<dim,spacedim>::face_iterator &face) const;
/**
- * Return equally spaced
- * intermediate points on a line.
+ * Return intermediate points on a line spaced according to the interior
+ * support points of the 1D Gauss-Lobatto quadrature formula.
*
- * The number of points requested
- * is given by the size of the
- * vector @p points. It is the
- * task of the derived classes to
- * arrange the points in
+ * The number of points requested is given by the size of the vector @p
+ * points. It is the task of the derived classes to arrange the points in
* approximately equal distances.
*
- * This function is called by the
- * @p MappingQ class. This
- * happens on each face line of a
- * cells that has got at least
- * one boundary line.
+ * This function is called by the @p MappingQ class. This happens on each
+ * face line of a cells that has got at least one boundary line.
*
- * As this function is not needed
- * for @p MappingQ1, it is not
- * made pure virtual, to avoid
- * the need to overload it. The
- * default implementation throws
- * an error in any case, however.
+ * As this function is not needed for @p MappingQ1, it is not made pure
+ * virtual, to avoid the need to overload it. The default implementation
+ * throws an error in any case, however.
*/
virtual
void
std::vector<Point<spacedim> > &points) const;
/**
- * Return equally spaced
- * intermediate points on a
- * boundary quad.
+ * Return intermediate points on a line spaced according to the tensor
+ * product of the interior support points of the 1D Gauss-Lobatto quadrature
+ * formula.
*
- * The number of points requested
- * is given by the size of the
- * vector @p points. It is
- * required that this number is a
- * square of another integer,
- * i.e. <tt>n=points.size()=m*m</tt>. It
- * is the task of the derived
- * classes to arrange the points
- * such they split the quad into
- * <tt>(m+1)(m+1)</tt> approximately
- * equal-sized subquads.
+ * The number of points requested is given by the size of the vector @p
+ * points. It is required that this number is a square of another integer,
+ * i.e. <tt>n=points.size()=m*m</tt>. It is the task of the derived classes
+ * to arrange the points such they split the quad into <tt>(m+1)(m+1)</tt>
+ * approximately equal-sized subquads.
*
- * This function is called by the
- * <tt>MappingQ<3></tt> class. This
- * happens each face quad of
- * cells in 3d that has got at
- * least one boundary face quad.
+ * This function is called by the <tt>MappingQ<3></tt> class. This happens
+ * each face quad of cells in 3d that has got at least one boundary face
+ * quad.
*
- * As this function is not needed
- * for @p MappingQ1, it is not
- * made pure virtual, to avoid
- * the need to overload it. The
- * default implementation throws
- * an error in any case, however.
+ * As this function is not needed for @p MappingQ1, it is not made pure
+ * virtual, to avoid the need to overload it. The default implementation
+ * throws an error in any case, however.
*/
virtual
void
std::vector<Point<spacedim> > &points) const;
/**
- * Depending on <tt>dim=2</tt> or
- * <tt>dim=3</tt> this function
- * calls the
- * get_intermediate_points_on_line
- * or the
- * get_intermediate_points_on_quad
- * function. It throws an
- * exception for
- * <tt>dim=1</tt>. This wrapper
- * allows dimension independent
- * programming.
+ * Depending on <tt>dim=2</tt> or <tt>dim=3</tt> this function calls the
+ * get_intermediate_points_on_line or the get_intermediate_points_on_quad
+ * function. It throws an exception for <tt>dim=1</tt>. This wrapper allows
+ * dimension independent programming.
*/
void
get_intermediate_points_on_face (const typename Triangulation<dim,spacedim>::face_iterator &face,
std::vector<Point<spacedim> > &points) const;
/**
- * Return the normal vector to the surface
- * at the point p. If p is not in fact
- * on the surface, but only closeby,
- * try to return something reasonable,
- * for example the normal vector
- * at the surface point closest to p.
- * (The point p will in fact not normally
- * lie on the actual surface, but rather
- * be a quadrature point mapped by some
- * polynomial mapping; the mapped surface,
- * however, will not usually coincide with
- * the actual surface.)
+ * Return the normal vector to the surface at the point p. If p is not in
+ * fact on the surface, but only close-by, try to return something
+ * reasonable, for example the normal vector at the surface point closest to
+ * p. (The point p will in fact not normally lie on the actual surface, but
+ * rather be a quadrature point mapped by some polynomial mapping; the
+ * mapped surface, however, will not usually coincide with the actual
+ * surface.)
*
- * The face iterator gives an indication
- * which face this function is supposed
- * to compute the normal vector for.
- * This is useful if the boundary of
- * the domain is composed of different
- * nondifferential pieces (for example
- * when using the StraightBoundary class
- * to approximate a geometry that is
- * completely described by the coarse mesh,
- * with piecewise (bi-)linear components
- * between the vertices, but where the
- * boundary may have a kink at the vertices
- * itself).
+ * The face iterator gives an indication which face this function is
+ * supposed to compute the normal vector for. This is useful if the
+ * boundary of the domain is composed of different nondifferential pieces
+ * (for example when using the StraightBoundary class to approximate a
+ * geometry that is completely described by the coarse mesh, with piecewise
+ * (bi-)linear components between the vertices, but where the boundary may
+ * have a kink at the vertices itself).
*
- * @note Implementations of this function
- * should be able to assume that the point p
- * lies within or close to the face described by the
- * first argument. In turn, callers of this
- * function should ensure that this is
- * in fact the case.
+ * @note Implementations of this function should be able to assume that the
+ * point p lies within or close to the face described by the first
+ * argument. In turn, callers of this function should ensure that this is in
+ * fact the case.
*/
virtual
Tensor<1,spacedim>
const Point<spacedim> &p) const;
/**
- * Compute the normal vectors to
- * the boundary at each vertex of
- * the given face. It is not
- * required that the normal
- * vectors be normed
- * somehow. Neither is it
- * required that the normals
- * actually point outward.
+ * Compute the normal vectors to the boundary at each vertex of the given
+ * face. It is not required that the normal vectors be normed
+ * somehow. Neither is it required that the normals actually point outward.
*
- * This function is
- * needed to compute data for C1
- * mappings. The default
- * implementation is to throw an
- * error, so you need not
- * overload this function in case
- * you do not intend to use C1
- * mappings.
+ * This function is needed to compute data for C1 mappings. The default
+ * implementation is to throw an error, so you need not overload this
+ * function in case you do not intend to use C1 mappings.
*
- * Note that when computing
- * normal vectors at a vertex
- * where the boundary is not
- * differentiable, you have to
- * make sure that you compute the
- * one-sided limits, i.e. limit
- * with respect to points inside
- * the given face.
+ * Note that when computing normal vectors at a vertex where the boundary is
+ * not differentiable, you have to make sure that you compute the one-sided
+ * limits, i.e. limit with respect to points inside the given face.
*/
virtual
void
FaceVertexNormals &face_vertex_normals) const;
/**
- * Given a candidate point and a
- * line segment characterized by
- * the iterator, return a point
- * that lies on the surface
- * described by this object. This
- * function is used in some mesh
- * smoothing algorithms that try
- * to move around points in order
- * to improve the mesh quality
- * but need to ensure that points
- * that were on the boundary
- * remain on the boundary.
+ * Given a candidate point and a line segment characterized by the iterator,
+ * return a point that lies on the surface described by this object. This
+ * function is used in some mesh smoothing algorithms that try to move
+ * around points in order to improve the mesh quality but need to ensure
+ * that points that were on the boundary remain on the boundary.
*
- * If spacedim==1, then the line
- * represented by the line
- * iterator is the entire space
- * (i.e. it is a cell, not a part
- * of the boundary), and the
- * returned point equals the
- * given input point.
+ * If spacedim==1, then the line represented by the line iterator is the
+ * entire space (i.e. it is a cell, not a part of the boundary), and the
+ * returned point equals the given input point.
*
- * Derived classes do not need to
- * implement this function unless
- * mesh smoothing algorithms are
- * used with a particular
- * boundary object. The default
- * implementation of this
- * function throws an exception
- * of type ExcPureFunctionCalled.
+ * Derived classes do not need to implement this function unless mesh
+ * smoothing algorithms are used with a particular boundary object. The
+ * default implementation of this function throws an exception of type
+ * ExcPureFunctionCalled.
*/
virtual
Point<spacedim>
const Point<spacedim> &candidate) const;
/**
- * Same function as above but for
- * a point that is to be
- * projected onto the area
- * characterized by the given
- * quad.
+ * Same function as above but for a point that is to be projected onto the
+ * area characterized by the given quad.
*
- * If spacedim<=2, then the surface
- * represented by the quad
- * iterator is the entire space
- * (i.e. it is a cell, not a part
- * of the boundary), and the
- * returned point equals the
- * given input point.
+ * If spacedim<=2, then the surface represented by the quad iterator is the
+ * entire space (i.e. it is a cell, not a part of the boundary), and the
+ * returned point equals the given input point.
*/
virtual
Point<spacedim>
const Point<spacedim> &candidate) const;
/**
- * Same function as above but for
- * a point that is to be
- * projected onto the area
- * characterized by the given
- * quad.
+ * Same function as above but for a point that is to be projected onto the
+ * area characterized by the given quad.
*
- * If spacedim<=3, then the manifold
- * represented by the hex
- * iterator is the entire space
- * (i.e. it is a cell, not a part
- * of the boundary), and the
- * returned point equals the
- * given input point.
+ * If spacedim<=3, then the manifold represented by the hex iterator is the
+ * entire space (i.e. it is a cell, not a part of the boundary), and the
+ * returned point equals the given input point.
*/
virtual
Point<spacedim>
project_to_surface (const typename Triangulation<dim,spacedim>::hex_iterator &hex,
const Point<spacedim> &candidate) const;
+
+protected:
+ /**
+ * Returns the support points of the Gauss-Lobatto quadrature formula used
+ * for intermediate points.
+ *
+ * @note Since the boundary description is closely tied to the unit cell
+ * support points of MappingQ, new boundary descriptions need to explicitly
+ * use these Gauss-Lobatto points and not equidistant points.
+ */
+ const std::vector<Point<1> > &
+ get_line_support_points (const unsigned int n_intermediate_points) const;
+
+private:
+ /**
+ * Point generator for the intermediate points on a boundary.
+ */
+ mutable std::vector<std_cxx1x::shared_ptr<QGaussLobatto<1> > > points;
+
+ /**
+ * Mutex for protecting the points array.
+ */
+ mutable Threads::Mutex mutex;
};
{
public:
/**
- * Default constructor. Some
- * compilers require this for
- * some reasons.
+ * Default constructor. Some compilers require this for some reasons.
*/
StraightBoundary ();
/**
- * Let the new point be the
- * arithmetic mean of the two
- * vertices of the line.
+ * Let the new point be the arithmetic mean of the two vertices of the line.
*
- * Refer to the general
- * documentation of this class
- * and the documentation of the
- * base class for more
- * information.
+ * Refer to the general documentation of this class and the documentation of
+ * the base class for more information.
*/
virtual Point<spacedim>
get_new_point_on_line (const typename Triangulation<dim,spacedim>::line_iterator &line) const;
/**
- * Let the new point be the
- * arithmetic mean of the four
- * vertices of this quad and the
- * four midpoints of the lines,
- * which are already created at
- * the time of calling this
- * function.
+ * Let the new point be the arithmetic mean of the four vertices of this
+ * quad and the four midpoints of the lines, which are already created at
+ * the time of calling this function.
*
- * Refer to the general
- * documentation of this class
- * and the documentation of the
- * base class for more
- * information.
+ * Refer to the general documentation of this class and the documentation of
+ * the base class for more information.
*/
virtual
Point<spacedim>
get_new_point_on_quad (const typename Triangulation<dim,spacedim>::quad_iterator &quad) const;
/**
- * Gives <tt>n=points.size()</tt>
- * points that splits the
- * StraightBoundary line into
- * $n+1$ partitions of equal
- * lengths.
+ * Gives <tt>n=points.size()</tt> points that splits the StraightBoundary
+ * line into $n+1$ partitions of equal lengths.
*
- * Refer to the general
- * documentation of this class
- * and the documentation of the
- * base class.
+ * Refer to the general documentation of this class and the documentation of
+ * the base class.
*/
virtual
void
std::vector<Point<spacedim> > &points) const;
/**
- * Gives <tt>n=points.size()=m*m</tt>
- * points that splits the
- * StraightBoundary quad into
- * $(m+1)(m+1)$ subquads of equal
- * size.
+ * Gives <tt>n=points.size()=m*m</tt> points that splits the
+ * StraightBoundary quad into $(m+1)(m+1)$ subquads of equal size.
*
- * Refer to the general
- * documentation of this class
- * and the documentation of the
- * base class.
+ * Refer to the general documentation of this class and the documentation of
+ * the base class.
*/
virtual
void
std::vector<Point<spacedim> > &points) const;
/**
- * Implementation of the function
- * declared in the base class.
+ * Implementation of the function declared in the base class.
*
- * Refer to the general
- * documentation of this class
- * and the documentation of the
- * base class.
+ * Refer to the general documentation of this class and the documentation of
+ * the base class.
*/
virtual
Tensor<1,spacedim>
const Point<spacedim> &p) const;
/**
- * Compute the normals to the
- * boundary at the vertices of
- * the given face.
+ * Compute the normals to the boundary at the vertices of the given face.
*
- * Refer to the general
- * documentation of this class
- * and the documentation of the
- * base class.
+ * Refer to the general documentation of this class and the documentation of
+ * the base class.
*/
virtual
void
typename Boundary<dim,spacedim>::FaceVertexNormals &face_vertex_normals) const;
/**
- * Given a candidate point and a
- * line segment characterized by
- * the iterator, return a point
- * that lies on the surface
- * described by this object. This
- * function is used in some mesh
- * smoothing algorithms that try
- * to move around points in order
- * to improve the mesh quality
- * but need to ensure that points
- * that were on the boundary
- * remain on the boundary.
+ * Given a candidate point and a line segment characterized by the iterator,
+ * return a point that lies on the surface described by this object. This
+ * function is used in some mesh smoothing algorithms that try to move
+ * around points in order to improve the mesh quality but need to ensure
+ * that points that were on the boundary remain on the boundary.
*
- * The point returned is the
- * projection of the candidate
- * point onto the line through
- * the two vertices of the given
- * line iterator.
+ * The point returned is the projection of the candidate point onto the line
+ * through the two vertices of the given line iterator.
*
- * If spacedim==1, then the line
- * represented by the line
- * iterator is the entire space
- * (i.e. it is a cell, not a part
- * of the boundary), and the
- * returned point equals the
- * given input point.
+ * If spacedim==1, then the line represented by the line iterator is the
+ * entire space (i.e. it is a cell, not a part of the boundary), and the
+ * returned point equals the given input point.
*/
virtual
Point<spacedim>
const Point<spacedim> &candidate) const;
/**
- * Same function as above but for
- * a point that is to be
- * projected onto the area
- * characterized by the given
- * quad.
+ * Same function as above but for a point that is to be projected onto the
+ * area characterized by the given quad.
*
- * The point returned is the
- * projection of the candidate
- * point onto the bilinear
- * surface spanned by the four
- * vertices of the given quad
- * iterator.
+ * The point returned is the projection of the candidate point onto the
+ * bilinear surface spanned by the four vertices of the given quad iterator.
*
- * If spacedim<=2, then the surface
- * represented by the quad
- * iterator is the entire space
- * (i.e. it is a cell, not a part
- * of the boundary), and the
- * returned point equals the
- * given input point.
+ * If spacedim<=2, then the surface represented by the quad iterator is the
+ * entire space (i.e. it is a cell, not a part of the boundary), and the
+ * returned point equals the given input point.
*/
virtual
Point<spacedim>
const Point<spacedim> &candidate) const;
/**
- * Same function as above but for
- * a point that is to be
- * projected onto the area
- * characterized by the given
- * quad.
+ * Same function as above but for a point that is to be projected onto the
+ * area characterized by the given quad.
*
- * The point returned is the
- * projection of the candidate
- * point onto the trilinear
- * manifold spanned by the eight
- * vertices of the given hex
+ * The point returned is the projection of the candidate point onto the
+ * trilinear manifold spanned by the eight vertices of the given hex
* iterator.
*
- * If spacedim<=3, then the manifold
- * represented by the hex
- * iterator is the entire space
- * (i.e. it is a cell, not a part
- * of the boundary), and the
- * returned point equals the
- * given input point.
+ * If spacedim<=3, then the manifold represented by the hex iterator is the
+ * entire space (i.e. it is a cell, not a part of the boundary), and the
+ * returned point equals the given input point.
*/
virtual
Point<spacedim>
const unsigned int dim = 2;
std::vector<Point<dim> > line_points (2);
- // 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. if not at the boundary,
- // get the respective points
- // from another function
+ // 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. if not at
+ // the boundary, get the respective points from another function
for (unsigned int line_no=0; line_no<GeometryInfo<dim>::lines_per_cell; ++line_no)
{
const Triangulation<dim>::line_iterator line = cell->line(line_no);
if (line->at_boundary())
{
- // first get the normal
- // vectors at the two
- // vertices of this line
- // from the boundary
- // description
+ // first get the normal vectors at the two vertices of this line
+ // from the boundary description
const Boundary<dim> &boundary
= line->get_triangulation().get_boundary(line->boundary_indicator());
Boundary<dim>::FaceVertexNormals face_vertex_normals;
boundary.get_normals_at_vertices (line, face_vertex_normals);
- // then transform them into
- // interpolation points for
- // a cubic polynomial
+ // then transform them into interpolation points for a cubic
+ // polynomial
//
- // for this, note that if
- // we describe the boundary
- // curve as a polynomial in
- // tangential coordinate
- // @p{t=0..1} (along the
- // line) and @p{s} in
- // normal direction, then
- // the cubic mapping is
- // such that @p{s = a*t**3
- // + b*t**2 + c*t + d}, and
- // we want to determine the
- // interpolation points at
- // @p{t=1/3} and
- // @p{t=2/3}. Since at
- // @p{t=0,1} we want a
- // vertex which is actually
- // at the boundary, we know
- // that @p{d=0} and
- // @p{a=-b-c}. As
- // side-conditions, we want
- // that the derivatives at
- // @p{t=0} and @p{t=1},
- // i.e. at the vertices
- // match those returned by
- // the boundary. We then
- // have that
- // @p{s(1/3)=1/27(2b+8c)}
- // and
- // @p{s(2/3)=4/27b+10/27c}.
+ // for this, note that if we describe the boundary curve as a
+ // polynomial in tangential coordinate @p{t=0..1} (along the line)
+ // and @p{s} in normal direction, then the cubic mapping is such
+ // that @p{s = a*t**3 + b*t**2 + c*t + d}, and we want to determine
+ // the interpolation points at @p{t=0.276} and @p{t=0.724}
+ // (Gauss-Lobatto points). Since at @p{t=0,1} we want a vertex which
+ // is actually at the boundary, we know that @p{d=0} and @p{a=-b-c},
+ // which gives @p{s(0.276)} and @{s(0.726)} in terms of @p{b,c}. As
+ // side-conditions, we want that the derivatives at @p{t=0} and
+ // @p{t=1}, i.e. at the vertices match those returned by the
+ // boundary.
//
- // The task is then first
- // to determine the
- // coefficients from the
- // tangentials. for that,
- // first rotate the
- // tangents of @p{s(t)}
- // into the global
- // coordinate system. they
- // are @p{A (1,c)} and @p{A
- // (1,-b-2c)} with @p{A} the
- // rotation matrix, since
- // the tangentials in the
- // coordinate system
- // relative to the line are
- // @p{(1,c)} and @p{(1,-b-2c)}
- // at the two vertices,
- // respectively. We then
- // have to make sure by
- // matching @p{b,c} that
- // these tangentials are
- // orthogonal to the normals
- // returned by the boundary
- // object
+ // The task is then first to determine the coefficients from the
+ // tangentials. for that, first rotate the tangents of @p{s(t)} into
+ // the global coordinate system. they are @p{A (1,c)} and @p{A
+ // (1,-b-2c)} with @p{A} the rotation matrix, since the tangentials
+ // in the coordinate system relative to the line are @p{(1,c)} and
+ // @p{(1,-b-2c)} at the two vertices, respectively. We then have to
+ // make sure by matching @p{b,c} that these tangentials are
+ // orthogonal to the normals returned by the boundary object
const Tensor<1,2> coordinate_vector = line->vertex(1) - line->vertex(0);
const double h = std::sqrt(coordinate_vector * coordinate_vector);
Tensor<1,2> coordinate_axis = coordinate_vector;
-face_vertex_normals[1][0] * std::sin(alpha)))
-2*c;
+ QGaussLobatto<1> quad_points(4);
+ const double t1 = quad_points.point(1)[0];
+ const double t2 = quad_points.point(2)[0];
+ const double s_t1 = (((-b-c)*t1+b)*t1+c)*t1;
+ const double s_t2 = (((-b-c)*t2+b)*t2+c)*t2;
- // next evaluate the so
- // determined cubic
- // polynomial at the points
- // 1/3 and 2/3, first in
- // unit coordinates
- const Point<2> new_unit_points[2] = { Point<2>(1./3., 1./27.*(2*b+8*c)),
- Point<2>(2./3., 4./27.*b+10./27.*c)
+ // next evaluate the so determined cubic polynomial at the points
+ // 1/3 and 2/3, first in unit coordinates
+ const Point<2> new_unit_points[2] = { Point<2>(t1, s_t1),
+ Point<2>(t2, s_t2)
};
- // then transform these
- // points to real
- // coordinates by rotating,
+ // then transform these points to real coordinates by rotating,
// scaling and shifting
for (unsigned int i=0; i<2; ++i)
{
// Qp mapping of cells at the boundary.
const QGaussLobatto<1> points(degree+1);
tensor_pols = new TensorProductPolynomials<dim>
- (Polynomials::LagrangeEquidistant::generate_complete_basis(degree));
+ (Polynomials::generate_complete_Lagrange_basis(points.get_points()));
Assert (n_shape_functions==tensor_pols->n(),
ExcInternalError());
Assert(n_inner+n_outer==n_shape_functions, ExcInternalError());
}
-// Note that the CellSimilarity flag is
-// modifyable, since MappingQ can need to
-// recalculate data even when cells are
-// similar.
+// Note that the CellSimilarity flag is modifyable, since MappingQ can need to
+// recalculate data even when cells are similar.
template<int dim, int spacedim>
void
MappingQ<dim,spacedim>::fill_fe_values (
std::vector<Point<spacedim> > &normal_vectors,
CellSimilarity::Similarity &cell_similarity) const
{
- // convert data object to internal
- // data for this class. fails with
- // an exception if that is not
- // possible
+ // convert data object to internal data for this class. fails with an
+ // exception if that is not possible
Assert (dynamic_cast<InternalData *> (&mapping_data) != 0, ExcInternalError());
InternalData &data = static_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
+ // 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. furthermore, we need to
- // ensure that the flag indicating
- // whether we can use some similarity has
- // to be modified - for a general
- // MappingQ, the data needs to be
- // recomputed anyway since then the
- // mapping changes the data. this needs
- // to be known also for later operations,
- // so modify the variable here. this also
- // affects the calculation of the next
- // cell -- if we use Q1 data on the next
- // cell, the data will still be invalid.
+ // depending on this result, use this or the other data object for the
+ // mapping. furthermore, we need to ensure that the flag indicating whether
+ // we can use some similarity has to be modified - for a general MappingQ,
+ // the data needs to be recomputed anyway since then the mapping changes the
+ // data. this needs to be known also for later operations, so modify the
+ // variable here. this also affects the calculation of the next cell -- if
+ // we use Q1 data on the next cell, the data will still be invalid.
typename MappingQ1<dim,spacedim>::InternalData *p_data=0;
if (data.use_mapping_q1_on_current_cell)
p_data=&data.mapping_q1_data;
std::vector<Tensor<1,spacedim> > &exterior_forms,
std::vector<Point<spacedim> > &normal_vectors) const
{
- // convert data object to internal
- // data for this class. fails with
- // an exception if that is not
- // possible
+ // convert data object to internal data for this class. fails with an
+ // exception if that is not possible
Assert (dynamic_cast<InternalData *> (&mapping_data) != 0,
ExcInternalError());
InternalData &data = static_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 entirely in
- // the interior of the domain. note
- // that it is not sufficient to ask
- // whether the present _face_ is in
- // the interior, as the mapping on
- // the face depends on the mapping
- // of the cell, which in turn
- // depends on the fact whether
- // _any_ of the faces of this cell
- // is at the boundary, not only the
- // present face
+ // check whether this cell needs the full mapping or can be treated by a
+ // reduced Q1 mapping, e.g. if the cell is entirely in the interior of the
+ // domain. note that it is not sufficient to ask whether the present _face_
+ // is in the interior, as the mapping on the face depends on the mapping of
+ // the cell, which in turn depends on the fact whether _any_ of the faces of
+ // this cell is at the boundary, not only the present face
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
+ // depending on this result, use this or the other data object for the
+ // mapping
typename MappingQ1<dim,spacedim>::InternalData *p_data=0;
if (data.use_mapping_q1_on_current_cell)
p_data=&data.mapping_q1_data;
std::vector<Tensor<1,spacedim> > &exterior_forms,
std::vector<Point<spacedim> > &normal_vectors) const
{
- // convert data object to internal
- // data for this class. fails with
- // an exception if that is not
- // possible
+ // convert data object to internal data for this class. fails with an
+ // exception if that is not possible
Assert (dynamic_cast<InternalData *> (&mapping_data) != 0,
ExcInternalError());
InternalData &data = static_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 entirely in
- // the interior of the domain. note
- // that it is not sufficient to ask
- // whether the present _face_ is in
- // the interior, as the mapping on
- // the face depends on the mapping
- // of the cell, which in turn
- // depends on the fact whether
- // _any_ of the faces of this cell
- // is at the boundary, not only the
- // present face
+ // check whether this cell needs the full mapping or can be treated by a
+ // reduced Q1 mapping, e.g. if the cell is entirely in the interior of the
+ // domain. note that it is not sufficient to ask whether the present _face_
+ // is in the interior, as the mapping on the face depends on the mapping of
+ // the cell, which in turn depends on the fact whether _any_ of the faces of
+ // this cell is at the boundary, not only the present face
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
+ // depending on this result, use this or the other data object for the
+ // mapping
typename MappingQ1<dim,spacedim>::InternalData *p_data=0;
if (data.use_mapping_q1_on_current_cell)
p_data=&data.mapping_q1_data;
const unsigned int n_inner_2d=(degree-1)*(degree-1);
const unsigned int n_outer_2d=4+4*(degree-1);
- // first check whether we have precomputed
- // the values for some polynomial degree;
- // the sizes of arrays is
- // n_inner_2d*n_outer_2d
+ // first check whether we have precomputed the values for some polynomial
+ // degree; the sizes of arrays is n_inner_2d*n_outer_2d
double const *loqv_ptr=0;
switch (degree)
{
- // for degree==1, we shouldn't have to
- // compute any support points, since
+ // for degree==1, we shouldn't have to compute any support points, since
// all of them are on the vertices
case 2:
{
- // (checked these values against the
- // output of compute_laplace_vector
- // again, and found they're indeed
- // right -- just in case someone
- // wonders where they come from --
- // WB)
+ // (checked these values against the output of compute_laplace_vector
+ // again, and found they're indeed right -- just in case someone wonders
+ // where they come from -- WB)
static const double loqv2[1*8]
= {1/16., 1/16., 1/16., 1/16., 3/16., 3/16., 3/16., 3/16.};
loqv_ptr=&loqv2[0];
break;
}
- case 3:
- {
- // (same as above)
- static const double loqv3[4*12]
- = {80/1053., 1/81., 1/81., 11/1053., 25/117., 44/351.,
- 7/117., 16/351., 25/117., 44/351., 7/117., 16/351.,
- 1/81., 80/1053., 11/1053., 1/81., 7/117., 16/351.,
- 25/117., 44/351., 44/351., 25/117., 16/351., 7/117.,
- 1/81., 11/1053., 80/1053., 1/81., 44/351., 25/117.,
- 16/351., 7/117., 7/117., 16/351., 25/117., 44/351.,
- 11/1053., 1/81., 1/81., 80/1053., 16/351., 7/117.,
- 44/351., 25/117., 16/351., 7/117., 44/351., 25/117.
- };
- Assert (sizeof(loqv3)/sizeof(loqv3[0]) ==
- n_inner_2d * n_outer_2d,
- ExcInternalError());
-
- loqv_ptr=&loqv3[0];
-
- break;
- }
-
-
- case 4:
- {
- static const double loqv4[9*16]
- = {0.07405921850311571, -0.001075744628905992,
- -0.001075744628906007, 0.001914292239071463,
- 0.2231273865431892, 0.1346851306015187,
- 0.03812914216116724, 0.02913160002633252,
- 0.02200737428129396, 0.01600835564431224,
- 0.2231273865431891, 0.1346851306015187,
- 0.03812914216116723, 0.02913160002633253,
- 0.02200737428129391, 0.01600835564431222,
-
- 0.00664803151334206, 0.006648031513342719,
- 0.002873452861657458, 0.002873452861657626,
- 0.07903572682584378, 0.05969238281250031,
- 0.03619864817415824, 0.07903572682584187,
- 0.0596923828124999, 0.03619864817415815,
- 0.1527716818820237, 0.2348152760709273,
- 0.152771681882024, 0.02496269311797778,
- 0.04081948955407129, 0.02496269311797789,
-
- -0.001075744628906923, 0.07405921850311589,
- 0.001914292239071339, -0.001075744628905884,
- 0.02913160002633509, 0.02200737428129395,
- 0.01600835564431229, 0.2231273865431878,
- 0.1346851306015183, 0.0381291421611672,
- 0.03812914216116729, 0.1346851306015185,
- 0.2231273865431898, 0.01600835564431217,
- 0.02200737428129394, 0.02913160002633262,
-
- 0.006648031513342073, 0.002873452861657473,
- 0.006648031513342726, 0.002873452861657636,
- 0.1527716818820238, 0.2348152760709273,
- 0.152771681882024, 0.02496269311797779,
- 0.04081948955407131, 0.0249626931179779,
- 0.07903572682584376, 0.05969238281250026,
- 0.03619864817415824, 0.07903572682584187,
- 0.0596923828124998, 0.0361986481741581,
-
- 0.01106770833333302, 0.01106770833333336,
- 0.01106770833333337, 0.01106770833333374,
- 0.06770833333333424, 0.1035156250000011,
- 0.0677083333333344, 0.06770833333333376,
- 0.103515624999999, 0.06770833333333399,
- 0.06770833333333422, 0.1035156250000009,
- 0.06770833333333436, 0.0677083333333337,
- 0.1035156249999988, 0.0677083333333339,
-
- 0.002873452861657185, 0.006648031513342362,
- 0.002873452861657334, 0.006648031513343038,
- 0.02496269311797779, 0.04081948955407401,
- 0.02496269311797788, 0.1527716818820234,
- 0.234815276070926, 0.1527716818820237,
- 0.03619864817415819, 0.05969238281250028,
- 0.07903572682584407, 0.03619864817415804,
- 0.05969238281249986, 0.0790357268258422,
-
- -0.001075744628906913, 0.00191429223907134,
- 0.07405921850311592, -0.001075744628905865,
- 0.03812914216116729, 0.1346851306015185,
- 0.2231273865431899, 0.01600835564431217,
- 0.02200737428129396, 0.02913160002633264,
- 0.02913160002633509, 0.02200737428129391,
- 0.01600835564431228, 0.2231273865431878,
- 0.1346851306015183, 0.03812914216116718,
-
- 0.002873452861657176, 0.002873452861657321,
- 0.006648031513342374, 0.006648031513343037,
- 0.03619864817415817, 0.05969238281250032,
- 0.07903572682584409, 0.03619864817415805,
- 0.05969238281249992, 0.07903572682584221,
- 0.02496269311797776, 0.04081948955407392,
- 0.02496269311797785, 0.1527716818820233,
- 0.2348152760709258, 0.1527716818820236,
-
- 0.001914292239071237, -0.001075744628906803,
- -0.001075744628906778, 0.07405921850311617,
- 0.01600835564431228, 0.02200737428129401,
- 0.02913160002633524, 0.03812914216116726,
- 0.1346851306015182, 0.2231273865431886,
- 0.01600835564431228, 0.02200737428129397,
- 0.02913160002633523, 0.03812914216116726,
- 0.1346851306015181, 0.2231273865431886,
- };
-
- Assert (sizeof(loqv4)/sizeof(loqv4[0]) ==
- n_inner_2d * n_outer_2d,
- ExcInternalError());
-
- loqv_ptr=&loqv4[0];
-
- break;
- }
-
- // no other cases implemented,
- // so simply fall through
+ // no other cases implemented, so simply fall through
default:
break;
}
if (loqv_ptr!=0)
{
- // precomputed. copy values to
- // the loqvs array
+ // precomputed. copy values to the loqvs array
loqvs.reinit(n_inner_2d, n_outer_2d);
for (unsigned int unit_point=0; unit_point<n_inner_2d; ++unit_point)
for (unsigned int k=0; k<n_outer_2d; ++k)
else
{
// not precomputed, then do so now
- if (dim==2)
+ if (dim == 2)
compute_laplace_vector(loqvs);
- else
- // computing the Laplace vector for
- // faces is not supported in 3d at
- // present. presumably, doing so
- // would not be so hard: we would
- // only have to call the function in
- // 2d, i.e. the quad(=face) values in
- // 3d are equal to the quad(=cell)
- // values in 2d. however, that would
- // require us to link in the 2d
- // library, which is kind of awkward
- // (note that compute_laplace_vector
- // really makes use of a lot of 2d
- // stuff, such as FEValues etc). an
- // alternative would be to precompute
- // the values of this array for a
- // couple of higher mapping orders,
- // pin down their values and insert
- // them into the array above.
- Assert (false, ExcNotImplemented());
+ else if (dim == 3)
+ {
+ MappingQ<2,2> mapping_2d(this->degree);
+ loqvs = mapping_2d.laplace_on_quad_vector;
+ }
}
- // the sum of weights of the points
- // at the outer rim should be
- // one. check this
+ // the sum of weights of the points at the outer rim should be one. check
+ // this
for (unsigned int unit_point=0; unit_point<loqvs.n_rows(); ++unit_point)
Assert(std::fabs(std::accumulate(loqvs[unit_point].begin(),
loqvs[unit_point].end(),0.)-1)<1e-13*this->degree,
{
Assert(degree>1, ExcInternalError());
- // first check whether we have
- // precomputed the values for some
- // polynomial degree
+ // first check whether we have precomputed the values for some polynomial
+ // degree
double const *lohv_ptr=0;
if (degree==2)
{
if (lohv_ptr!=0)
{
- // precomputed. copy values to
- // the lohvs array
+ // precomputed. copy values to the lohvs array
lohvs.reinit(n_inner, n_outer);
for (unsigned int unit_point=0; unit_point<n_inner; ++unit_point)
for (unsigned int k=0; k<n_outer; ++k)
// not precomputed, then do so now
compute_laplace_vector(lohvs);
- // the sum of weights of the points
- // at the outer rim should be
- // one. check this
+ // the sum of weights of the points at the outer rim should be one. check
+ // this
for (unsigned int unit_point=0; unit_point<n_inner; ++unit_point)
Assert(std::fabs(std::accumulate(lohvs[unit_point].begin(),
lohvs[unit_point].end(),0.) - 1)<1e-12*this->degree*this->degree,
Assert(lvs.n_rows()==0, ExcInternalError());
Assert(dim==2 || dim==3, ExcNotImplemented());
- // for degree==1, we shouldn't have to
- // compute any support points, since all of
- // them are on the vertices
+ // for degree==1, we shouldn't have to compute any support points, since all
+ // of them are on the vertices
Assert(degree>1, ExcInternalError());
- // compute the shape
- // gradients at the quadrature
- // points on the unit cell
+ // compute the shape gradients at the quadrature points on the unit cell
const QGauss<dim> quadrature(degree+1);
const unsigned int n_q_points=quadrature.size();
quadrature_data.shape_derivatives.resize(n_shape_functions * n_q_points);
this->compute_shapes(quadrature.get_points(), quadrature_data);
- // Compute the stiffness matrix of
- // the inner dofs
+ // Compute the stiffness matrix of the inner dofs
FullMatrix<long double> S(n_inner);
for (unsigned int point=0; point<n_q_points; ++point)
for (unsigned int i=0; i<n_inner; ++i)
S(i,j) += res * (long double)quadrature.weight(point);
}
- // Compute the components of T to be the
- // product of gradients of inner and
+ // Compute the components of T to be the product of gradients of inner and
// outer shape functions.
FullMatrix<long double> T(n_inner, n_outer);
for (unsigned int point=0; point<n_q_points; ++point)
// S:=S_1*T
S_1.mmult(S_1_T,T);
- // Resize and initialize the
- // lvs
+ // Resize and initialize the lvs
lvs.reinit (n_inner, n_outer);
for (unsigned int i=0; i<n_inner; ++i)
for (unsigned int k=0; k<n_outer; ++k)
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
+ // 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(degree));
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
+ // 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());
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 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
+ // 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,spacedim>::compute_support_points_laplace(const typename Triangulation<dim,spacedim>::cell_iterator &cell,
std::vector<Point<spacedim> > &a) const
{
- // in any case, we need the
- // vertices first
+ // 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);
add_line_support_points(cell, a);
break;
case 2:
- // in 2d, add the
- // points on the four
- // bounding lines to
- // the exterior (outer)
- // points
+ // 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);
case 3:
{
- // in 3d also add the
- // points located on
- // the boundary faces
+ // 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);
MappingQ<1>::add_line_support_points (const Triangulation<1>::cell_iterator &,
std::vector<Point<1> > &) const
{
- // there are no points on bounding
- // lines which are to be added
+ // there are no points on bounding lines which are to be added
const unsigned int dim=1;
Assert (dim > 1, ExcImpossibleInDim(dim));
}
{
const unsigned int dim = 1;
const unsigned int spacedim = 2;
- // Ask for the mid point, if that's
- // the only thing we need.
+ // Ask for the mid point, if that's the only thing we need.
if (degree==2)
{
const Boundary<dim,spacedim> *const boundary
a.push_back(boundary->get_new_point_on_line(cell));
}
else
- // otherwise call the more
- // complicated functions and ask
- // for inner points from the
- // boundary description
+ // otherwise call the more complicated functions and ask for inner points
+ // from the boundary description
{
std::vector<Point<spacedim> > line_points (degree-1);
{
const unsigned int dim = 1;
const unsigned int spacedim = 3;
- // Ask for the mid point, if that's
- // the only thing we need.
+ // Ask for the mid point, if that's the only thing we need.
if (degree==2)
{
const Boundary<dim,spacedim> *const boundary
a.push_back(boundary->get_new_point_on_line(cell));
}
else
- // otherwise call the more
- // complicated functions and ask
- // for inner points from the
- // boundary description
+ // otherwise call the more complicated functions and ask for inner points
+ // from the boundary description
{
std::vector<Point<spacedim> > line_points (degree-1);
std::vector<Point<spacedim> > &a) const
{
static const StraightBoundary<dim,spacedim> straight_boundary;
- // if we only need the midpoint,
- // then ask for it.
+ // if we only need the midpoint, then ask for it.
if (degree==2)
{
for (unsigned int line_no=0; line_no<GeometryInfo<dim>::lines_per_cell; ++line_no)
};
}
else
- // otherwise call the more
- // complicated functions and ask
- // for inner points from the
- // boundary description
+ // otherwise call the more complicated functions and ask for inner points
+ // from the boundary description
{
std::vector<Point<spacedim> > 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
+ // 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,spacedim>::line_iterator line = cell->line(line_no);
boundary->get_intermediate_points_on_line (line, line_points);
if (dim==3)
{
- // in 3D, lines might be in wrong
- // orientation. if so, reverse
+ // in 3D, lines might be in wrong orientation. if so, reverse
// the vector
if (cell->line_orientation(line_no))
a.insert (a.end(), line_points.begin(), line_points.end());
a.insert (a.end(), line_points.rbegin(), line_points.rend());
}
else
- // in 2D, lines always have the
- // correct orientation. simply
- // append all points
+ // in 2D, lines always have the correct orientation. simply append
+ // all points
a.insert (a.end(), line_points.begin(), line_points.end());
}
vertices_per_cell = GeometryInfo<3>::vertices_per_cell;
static const StraightBoundary<3> straight_boundary;
- // used if face quad at boundary or
- // entirely in the interior of the
- // domain
+ // used if face quad at boundary or entirely in the interior of the domain
std::vector<Point<3> > quad_points ((degree-1)*(degree-1));
- // used if only one line of face
- // quad is at boundary
+ // used if only one line of face quad is at boundary
std::vector<Point<3> > b(4*degree);
- // loop over all faces and collect
- // points on them
+ // loop over all faces and collect points on them
for (unsigned int face_no=0; face_no<faces_per_cell; ++face_no)
{
const Triangulation<3>::face_iterator face = cell->face(face_no);
- // select the correct mappings
- // for the present face
+ // select the correct mappings for the present face
const bool face_orientation = cell->face_orientation(face_no),
face_flip = cell->face_flip (face_no),
face_rotation = cell->face_rotation (face_no);
face_rotation)),
ExcInternalError());
- // indices of the lines that
- // bound a face are given by
- // GeometryInfo<3>::
+ // indices of the lines that bound a face are given by GeometryInfo<3>::
// face_to_cell_lines
for (unsigned int i=0; i<lines_per_face; ++i)
Assert(face->line(i)==cell->line(GeometryInfo<3>::face_to_cell_lines(
ExcInternalError());
#endif
- // if face at boundary, then
- // ask boundary object to
- // return intermediate points
- // on it
+ // if face at boundary, then ask boundary object to return intermediate
+ // points on it
if (face->at_boundary())
{
face->get_triangulation().get_boundary(face->boundary_indicator())
.get_intermediate_points_on_quad (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.
+ // 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.
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,
}
else
{
- // face is not at boundary,
- // but maybe some of its
- // lines are. count them
+ // face is not at boundary, but maybe some of its lines are. count
+ // them
unsigned int lines_at_boundary=0;
for (unsigned int i=0; i<lines_per_face; ++i)
if (face->line(i)->at_boundary())
Assert(lines_at_boundary<=lines_per_face, ExcInternalError());
- // if at least one of the
- // lines bounding this quad
- // is at the boundary, then
- // collect points
- // separately
+ // if at least one of the lines bounding this quad is at the
+ // boundary, then collect points separately
if (lines_at_boundary>0)
{
- // call of function
- // apply_laplace_vector
- // increases size of b
- // about 1. There
- // resize b for the
- // case the mentioned
- // function was already
- // called.
+ // call of function apply_laplace_vector increases size of b
+ // about 1. There resize b for the case the mentioned function
+ // was already called.
b.resize(4*degree);
- // b is of size
- // 4*degree, make sure
- // that this is the
- // right size
+ // b is of size 4*degree, make sure that this is the right size
Assert(b.size()==vertices_per_face+lines_per_face*(degree-1),
ExcDimensionMismatch(b.size(),
vertices_per_face+lines_per_face*(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
- // orientation. Doing so, the
- // calculated points will be in
+ // 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
+ // orientation. Doing so, the calculated points will be in
// standard orientation as well.
for (unsigned int i=0; i<vertices_per_face; ++i)
b[i]=a[GeometryInfo<3>::face_to_cell_vertices(face_no, i)];
a[vertices_per_cell + GeometryInfo<3>::face_to_cell_lines(
face_no, i)*(degree-1)+j];
- // Now b includes the support
- // points on the quad and we can
+ // 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*degree+(degree-1)*(degree-1),
}
else
{
- // face is entirely in
- // the interior. get
- // intermediate points
- // from a straight
- // boundary object
+ // face is entirely in the interior. get intermediate points
+ // from a straight boundary object
straight_boundary.get_intermediate_points_on_quad (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.
+ // 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.
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,
const MappingType mapping_type) const
{
AssertDimension (input.size(), output.size());
- // The data object may be just a
- // MappingQ1::InternalData, so we
- // have to test for this first.
+ // The data object may be just a MappingQ1::InternalData, so we have to test
+ // for this first.
const typename MappingQ1<dim,spacedim>::InternalData *q1_data =
dynamic_cast<const typename MappingQ1<dim,spacedim>::InternalData *> (&mapping_data);
Assert(q1_data!=0, ExcInternalError());
- // If it is a genuine
- // MappingQ::InternalData, we have
- // to test further
+ // If it is a genuine MappingQ::InternalData, we have to test further
if (!q1_data->is_mapping_q1_data)
{
Assert (dynamic_cast<const InternalData *>(&mapping_data) != 0,
ExcInternalError());
const InternalData &data = static_cast<const InternalData &>(mapping_data);
- // If we only use the
- // Q1-portion, we have to
- // extract that data object
+ // If we only use the Q1-portion, we have to extract that data object
if (data.use_mapping_q1_on_current_cell)
q1_data = &data.mapping_q1_data;
}
- // Now, q1_data should have the
- // right tensors in it and we call
- // the base classes transform
- // function
+ // Now, q1_data should have the right tensors in it and we call the base
+ // classes transform function
MappingQ1<dim,spacedim>::transform(input, output, *q1_data, mapping_type);
}
const MappingType mapping_type) const
{
AssertDimension (input.size(), output.size());
- // The data object may be just a
- // MappingQ1::InternalData, so we
- // have to test for this first.
+ // The data object may be just a MappingQ1::InternalData, so we have to test
+ // for this first.
const typename MappingQ1<dim,spacedim>::InternalData *q1_data =
dynamic_cast<const typename MappingQ1<dim,spacedim>::InternalData *> (&mapping_data);
Assert(q1_data!=0, ExcInternalError());
- // If it is a genuine
- // MappingQ::InternalData, we have
- // to test further
+ // If it is a genuine MappingQ::InternalData, we have to test further
if (!q1_data->is_mapping_q1_data)
{
Assert (dynamic_cast<const InternalData *>(&mapping_data) != 0,
ExcInternalError());
const InternalData &data = static_cast<const InternalData &>(mapping_data);
- // If we only use the
- // Q1-portion, we have to
- // extract that data object
+ // If we only use the Q1-portion, we have to extract that data object
if (data.use_mapping_q1_on_current_cell)
q1_data = &data.mapping_q1_data;
}
- // Now, q1_data should have the
- // right tensors in it and we call
- // the base classes transform
- // function
+ // Now, q1_data should have the right tensors in it and we call the base
+ // classes transform function
MappingQ1<dim,spacedim>::transform(input, output, *q1_data, mapping_type);
}
const MappingType mapping_type) const
{
AssertDimension (input.size(), output.size());
- // The data object may be just a
- // MappingQ1::InternalData, so we
- // have to test for this first.
+ // The data object may be just a MappingQ1::InternalData, so we have to test
+ // for this first.
const typename MappingQ1<dim,spacedim>::InternalData *q1_data =
dynamic_cast<const typename MappingQ1<dim,spacedim>::InternalData *> (&mapping_data);
Assert(q1_data!=0, ExcInternalError());
- // If it is a genuine
- // MappingQ::InternalData, we have
- // to test further
+ // If it is a genuine MappingQ::InternalData, we have to test further
if (!q1_data->is_mapping_q1_data)
{
Assert (dynamic_cast<const InternalData *>(&mapping_data) != 0,
ExcInternalError());
const InternalData &data = static_cast<const InternalData &>(mapping_data);
- // If we only use the
- // Q1-portion, we have to
- // extract that data object
+ // If we only use the Q1-portion, we have to extract that data object
if (data.use_mapping_q1_on_current_cell)
q1_data = &data.mapping_q1_data;
}
- // Now, q1_data should have the
- // right tensors in it and we call
- // the base classes transform
- // function
+ // Now, q1_data should have the right tensors in it and we call the base
+ // classes transform function
MappingQ1<dim,spacedim>::transform(input, output, *q1_data, mapping_type);
}
transform_unit_to_real_cell (const typename Triangulation<dim,spacedim>::cell_iterator &cell,
const Point<dim> &p) const
{
- // Use the get_data function to
- // create an InternalData with data
- // vectors of the right size and
- // transformation shape values
- // already computed at point p.
+ // Use the get_data function to create an InternalData with data vectors of
+ // the right size and transformation shape values already computed at point
+ // p.
const Quadrature<dim> point_quadrature(p);
std::auto_ptr<InternalData>
mdata (dynamic_cast<InternalData *> (
&*mdata);
compute_mapping_support_points(cell, p_data->mapping_support_points);
- // If this should be Q1, ignore all
- // other support points.
+ // If this should be Q1, ignore all other support points.
if (p_data->shape_values.size()<p_data->mapping_support_points.size())
p_data->mapping_support_points.resize
(GeometryInfo<dim>::vertices_per_cell);
transform_real_to_unit_cell (const typename Triangulation<dim,spacedim>::cell_iterator &cell,
const Point<spacedim> &p) const
{
- // first a Newton iteration based
- // on a Q1 mapping to get a good
- // starting point, the idea being
- // that this is cheaper than trying
- // to start with the real mapping
- // and likely also more robust.
+ // first a Newton iteration based on a Q1 mapping to get a good starting
+ // point, the idea being that this is cheaper than trying to start with the
+ // real mapping and likely also more robust.
//
- // that said, this doesn't always
- // work: there are cases where the
- // point is outside the cell and
- // the inverse mapping doesn't
- // converge. in that case, use the
- // center point of the cell as a
- // starting point if we are to go
- // on using the full mapping, or
- // just propagate up the exception
- // if we had no intention of
- // continuing with the full mapping
+ // that said, this doesn't always work: there are cases where the point is
+ // outside the cell and the inverse mapping doesn't converge. in that case,
+ // use the center point of the cell as a starting point if we are to go on
+ // using the full mapping, or just propagate up the exception if we had no
+ // intention of continuing with the full mapping
Point<dim> initial_p_unit;
try
{
}
catch (const typename Mapping<dim,spacedim>::ExcTransformationFailed &)
{
- // mirror the conditions of the
- // code below to determine if
- // we need to use an arbitrary
- // starting point or if we just
- // need to rethrow the
+ // mirror the conditions of the code below to determine if we need to
+ // use an arbitrary starting point or if we just need to rethrow the
// exception
if (cell->has_boundary_lines()
||
throw;
}
- // then a Newton iteration based on the
- // full MappingQ if we need this. note that
- // for interior cells with dim==spacedim,
- // the mapping used is in fact a Q1
- // mapping, so there is nothing we need to
- // do unless the iteration above failed
+ // then a Newton iteration based on the full MappingQ if we need this. note
+ // that for interior cells with dim==spacedim, the mapping used is in fact a
+ // Q1 mapping, so there is nothing we need to do unless the iteration above
+ // failed
if (cell->has_boundary_lines()
||
use_mapping_q_on_all_cells
||
(dim!=spacedim) )
{
- // use the full mapping. in case the
- // function above should have given us
- // something back that lies outside the
- // unit cell (that might happen because
- // we may have given a point 'p' that
- // lies inside the cell with the higher
- // order mapping, but outside the
- // Q1-mapped reference cell), then
- // project it back into the reference
- // cell in hopes that this gives a
- // better starting point to the
- // following iteration
+ // use the full mapping. in case the function above should have given us
+ // something back that lies outside the unit cell (that might happen
+ // because we may have given a point 'p' that lies inside the cell with
+ // the higher order mapping, but outside the Q1-mapped reference cell),
+ // then project it back into the reference cell in hopes that this gives
+ // a better starting point to the following iteration
initial_p_unit = GeometryInfo<dim>::project_to_unit_cell(initial_p_unit);
const Quadrature<dim> point_quadrature(initial_p_unit);
compute_mapping_support_points (cell, mdata->mapping_support_points);
- // If this is a q1 mapping,
- // then only use the support
- // points on the vertices.
+ // If this is a q1 mapping, then only use the support points on the
+ // vertices.
if (mdata->shape_values.size() < mdata->mapping_support_points.size())
mdata->mapping_support_points.resize(GeometryInfo<dim>::vertices_per_cell);
-
return this->transform_real_to_unit_cell_internal(cell, p, initial_p_unit, *mdata);
}
else
:
Quadrature<dim>(Utilities::fixed_power<dim>(map_degree+1))
{
- // first we determine the support points
- // on the unit cell in lexicographic order.
- // for this purpose we can use an interated
- // trapezoidal quadrature rule.
- const QTrapez<1> q1d;
- const QIterated<dim> q_iterated(q1d,map_degree);
+ // first we determine the support points on the unit cell in lexicographic
+ // order, which are (in accordance with MappingQ) the support points of
+ // QGaussLobatto.
+ const QGaussLobatto<dim> q_iterated(map_degree+1);
const unsigned int n_q_points = q_iterated.size();
- // we then need to define a renumbering
- // vector that allows us to go from a
- // lexicographic numbering scheme to a hierarchic
- // one. this fragment is taking almost verbatim
- // from the MappingQ class.
-
+ // we then need to define a renumbering vector that allows us to go from a
+ // lexicographic numbering scheme to a hierarchic one. this fragment is
+ // taking almost verbatim from the MappingQ class.
std::vector<unsigned int> renumber(n_q_points);
std::vector<unsigned int> dpo(dim+1, 1U);
for (unsigned int i=1; i<dpo.size(); ++i)
FETools::lexicographic_to_hierarchic_numbering (
FiniteElementData<dim> (dpo, 1, map_degree), renumber);
- // finally we assign the quadrature points in the
- // required order.
-
+ // finally we assign the quadrature points in the required order.
for (unsigned int q=0; q<n_q_points; ++q)
this->quadrature_points[renumber[q]] = q_iterated.point(q);
}
std::vector<Point<spacedim> > &a) const
{
- // first, basic assertion
- // with respect to vector size,
+ // first, basic assertion with respect to vector size,
const types::global_dof_index n_dofs = euler_dof_handler->n_dofs();
const types::global_dof_index vector_size = euler_vector->size();
AssertDimension(vector_size,n_dofs);
- // we then transform our tria iterator
- // into a dof iterator so we can
- // access data not associated with
- // triangulations
+ // we then transform our tria iterator into a dof iterator so we can access
+ // data not associated with triangulations
typename DoFHandler<dim,spacedim>::cell_iterator dof_cell(*cell, euler_dof_handler);
Assert (dof_cell->active() == true, ExcInactiveCell());
- // our quadrature rule is chosen
- // so that each quadrature point
- // corresponds to a support point
- // in the undeformed configuration.
- // we can then query the given
- // displacement field at these points
- // to determine the shift vector that
- // maps the support points to the
- // deformed configuration.
-
- // we assume that the given field contains
- // dim displacement components, but
- // that there may be other solution
- // components as well (e.g. pressures).
- // this class therefore assumes that the
- // first dim components represent the
- // actual shift vector we need, and simply
- // ignores any components after that.
- // this implies that the user should order
- // components appropriately, or create a
- // separate dof handler for the displacements.
+ // our quadrature rule is chosen so that each quadrature point corresponds
+ // to a support point in the undeformed configuration. we can then query
+ // the given displacement field at these points to determine the shift
+ // vector that maps the support points to the deformed configuration.
+
+ // we assume that the given field contains dim displacement components, but
+ // that there may be other solution components as well (e.g. pressures).
+ // this class therefore assumes that the first dim components represent the
+ // actual shift vector we need, and simply ignores any components after
+ // that. this implies that the user should order components appropriately,
+ // or create a separate dof handler for the displacements.
const unsigned int n_support_pts = support_quadrature.size();
const unsigned int n_components = euler_dof_handler->get_fe().n_components();
std::vector<Vector<double> > shift_vector(n_support_pts,Vector<double>(n_components));
- // fill shift vector for each
- // support point using an fe_values
- // object. make sure that the
- // fe_values variable isn't used
- // simulatenously from different
+ // fill shift vector for each support point using an fe_values object. make
+ // sure that the fe_values variable isn't used simulatenously from different
// threads
Threads::Mutex::ScopedLock lock(fe_values_mutex);
fe_values.reinit(dof_cell);
fe_values.get_function_values(*euler_vector, shift_vector);
- // and finally compute the positions of the
- // support points in the deformed
+ // and finally compute the positions of the support points in the deformed
// configuration.
-
a.resize(n_support_pts);
for (unsigned int q=0; q<n_support_pts; ++q)
{
std::vector<Point<spacedim> > &normal_vectors,
CellSimilarity::Similarity &cell_similarity) const
{
- // disable any previously detected
- // similarity and hand on to the respective
+ // disable any previously detected similarity and hand on to the respective
// function of the base class.
cell_similarity = CellSimilarity::invalid_next_cell;
MappingQ<dim,spacedim>::fill_fe_values (cell, q, mapping_data,
}
+
+template <int dim, int spacedim>
+const std::vector<Point<1> > &
+Boundary<dim,spacedim>::
+get_line_support_points (const unsigned int n_intermediate_points) const
+{
+ if (points.size() <= n_intermediate_points ||
+ points[n_intermediate_points].get() == 0)
+ {
+ Threads::Mutex::ScopedLock lock(mutex);
+ if (points.size() <= n_intermediate_points)
+ points.resize(n_intermediate_points+1);
+
+ // another thread might have created points in the meantime
+ if (points[n_intermediate_points].get() == 0)
+ {
+ std_cxx1x::shared_ptr<QGaussLobatto<1> >
+ quadrature (new QGaussLobatto<1>(n_intermediate_points+2));
+ points[n_intermediate_points] = quadrature;
+ }
+ }
+ return points[n_intermediate_points]->get_points();
+}
+
+
+
+
/* -------------------------- StraightBoundary --------------------- */
const unsigned int n=points.size();
Assert(n>0, ExcInternalError());
- const double dx=1./(n+1);
- double x=dx;
-
+ // Use interior points of QGaussLobatto quadrature formula support points
+ // for consistency with MappingQ
+ const std::vector<Point<1> > &line_points = this->get_line_support_points(n);
const Point<spacedim> vertices[2] = { line->vertex(0),
line->vertex(1)
};
- for (unsigned int i=0; i<n; ++i, x+=dx)
- points[i] = (1-x)*vertices[0] + x*vertices[1];
+ for (unsigned int i=0; i<n; ++i)
+ {
+ const double x = line_points[i+1][0];
+ points[i] = (1-x)*vertices[0] + x*vertices[1];
+ }
}
const unsigned int n=points.size();
Assert(n>0, ExcInternalError());
- const double dx=1./(n+1);
- double x=dx;
+ // Use interior points of QGaussLobatto quadrature formula support points
+ // for consistency with MappingQ
+ const std::vector<Point<1> > &line_points = this->get_line_support_points(n);
const Point<spacedim> vertices[2] = { line->vertex(0),
line->vertex(1)
};
- for (unsigned int i=0; i<n; ++i, x+=dx)
- points[i] = (1-x)*vertices[0] + x*vertices[1];
+ for (unsigned int i=0; i<n; ++i)
+ {
+ const double x = line_points[1+i][0];
+ points[i] = (1-x)*vertices[0] + x*vertices[1];
+ }
}
// is n a square number
Assert(m*m==n, ExcInternalError());
- const double ds=1./(m+1);
- double y=ds;
+ const std::vector<Point<1> > &line_points = this->get_line_support_points(m);
const Point<spacedim> vertices[4] = { quad->vertex(0),
quad->vertex(1),
quad->vertex(3)
};
- for (unsigned int i=0; i<m; ++i, y+=ds)
+ for (unsigned int i=0; i<m; ++i)
{
- double x=ds;
- for (unsigned int j=0; j<m; ++j, x+=ds)
- points[i*m+j]=((1-x) * vertices[0] +
- x * vertices[1]) * (1-y) +
- ((1-x) * vertices[2] +
- x * vertices[3]) * y;
+ const double y=line_points[1+i][0];
+ for (unsigned int j=0; j<m; ++j)
+ {
+ const double x=line_points[1+j][0];
+ points[i*m+j]=((1-x) * vertices[0] +
+ x * vertices[1]) * (1-y) +
+ ((1-x) * vertices[2] +
+ x * vertices[3]) * y;
+ }
}
}
// is n a square number
Assert(m*m==n, ExcInternalError());
- const double ds=1./(m+1);
- double y=ds;
+ const std::vector<Point<1> > &line_points = this->get_line_support_points(m);
const Point<spacedim> vertices[4] = { quad->vertex(0),
quad->vertex(1),
quad->vertex(3)
};
- for (unsigned int i=0; i<m; ++i, y+=ds)
+ for (unsigned int i=0; i<m; ++i)
{
- double x=ds;
- for (unsigned int j=0; j<m; ++j, x+=ds)
- points[i*m+j]=((1-x) * vertices[0] +
- x * vertices[1]) * (1-y) +
- ((1-x) * vertices[2] +
- x * vertices[3]) * y;
+ const double y=line_points[1+i][0];
+ for (unsigned int j=0; j<m; ++j)
+ {
+ const double x=line_points[1+j][0];
+ points[i*m+j]=((1-x) * vertices[0] +
+ x * vertices[1]) * (1-y) +
+ ((1-x) * vertices[2] +
+ x * vertices[3]) * y;
+ }
}
}
const unsigned int n=points.size();
Assert(n>0, ExcInternalError());
- // Do a simple linear interpolation
- // followed by projection, using
- // the same algorithm as above
- const Point<spacedim> ds = (v1-v0) / (n+1);
+ // Do a simple linear interpolation followed by projection, using the same
+ // algorithm as above
+ const std::vector<Point<1> > &line_points = this->get_line_support_points(n);
for (unsigned int i=0; i<n; ++i)
{
- const Point<spacedim> middle = v0 + (i+1)*ds;
+ const double x = line_points[i+1][0];
+ const Point<spacedim> middle = (1-x)*v0 + x*v1;
const Point<spacedim> vector_from_axis = (middle-point_on_axis) -
((middle-point_on_axis) * direction) * direction;
{
const unsigned int n = points.size ();
const Point<dim> axis = x_1 - x_0;
- const Point<dim> dx = (p1 - p0) / (n + 1);
Assert (n > 0, ExcInternalError ());
- for (unsigned int i = 0; i < n; ++i)
+ const std::vector<Point<1> > &line_points = this->get_line_support_points(n);
+
+ for (unsigned int i=0; i<n; ++i)
{
+ const double x = line_points[i+1][0];
+
// Compute the current point.
- const Point<dim> x_i = p0 + (i + 1) * dx;
- // To project this point on the
- // boundary of the cone we first
- // compute the orthogonal
- // projection of this point onto
- // the axis of the cone.
+ const Point<dim> x_i = (1-x)*p0 + x*p1;
+ // To project this point on the boundary of the cone we first compute
+ // the orthogonal projection of this point onto the axis of the cone.
const double c = (x_i - x_0) * axis / axis.square ();
const Point<dim> x_ip = x_0 + c * axis;
- // Compute the projection of
- // the middle point on the
- // boundary of the cone.
+ // Compute the projection of the middle point on the boundary of the
+ // cone.
points[i] = x_ip + get_radius (x_ip) * (x_i - x_ip) / (x_i - x_ip).norm ();
}
}
const Point<dim> axis = x_1 - x_0;
// Compute the middle point of the line.
const Point<dim> middle = StraightBoundary<dim>::get_new_point_on_line (line);
- // To project it on the boundary of
- // the cone we first compute the
- // orthogonal projection of the
- // middle point onto the axis of
- // the cone.
+ // To project it on the boundary of the cone we first compute the orthogonal
+ // projection of the middle point onto the axis of the cone.
const double c = (middle - x_0) * axis / axis.square ();
const Point<dim> middle_p = x_0 + c * axis;
- // Compute the projection of the
- // middle point on the boundary
- // of the cone.
+ // Compute the projection of the middle point on the boundary of the cone.
return middle_p + get_radius (middle_p) * (middle - middle_p) / (middle - middle_p).norm ();
}
const int dim = 3;
const Point<dim> axis = x_1 - x_0;
- // Compute the middle point of the
- // quad.
+ // Compute the middle point of the quad.
const Point<dim> middle = StraightBoundary<3>::get_new_point_on_quad (quad);
- // Same algorithm as above: To
- // project it on the boundary of
- // the cone we first compute the
- // orthogonal projection of the
- // middle point onto the axis of
- // the cone.
+ // Same algorithm as above: To project it on the boundary of the cone we
+ // first compute the orthogonal projection of the middle point onto the axis
+ // of the cone.
const double c = (middle - x_0) * axis / axis.square ();
const Point<dim> middle_p = x_0 + c * axis;
- // Compute the projection of the
- // middle point on the boundary
- // of the cone.
+ // Compute the projection of the middle point on the boundary of the cone.
return middle_p + get_radius (middle_p) * (middle - middle_p) / (middle - middle_p).norm ();
}
for (unsigned int vertex = 0; vertex < GeometryInfo<dim>::vertices_per_cell; ++vertex)
{
- // Compute the orthogonal
- // projection of the vertex onto
- // the axis of the cone.
+ // Compute the orthogonal projection of the vertex onto the axis of the
+ // cone.
const double c = (face->vertex (vertex) - x_0) * axis / axis.square ();
const Point<dim> vertex_p = x_0 + c * axis;
- // Then compute the vector
- // pointing from the point
- // <tt>vertex_p</tt> on the axis
- // to the vertex.
+ // Then compute the vector pointing from the point <tt>vertex_p</tt> on
+ // the axis to the vertex.
const Point<dim> axis_to_vertex = face->vertex (vertex) - vertex_p;
face_vertex_normals[vertex] = axis_to_vertex / axis_to_vertex.norm ();
Assert(std::fabs(v1.square()-r2)<eps*r2, ExcInternalError());
const double alpha=std::acos((v0*v1)/std::sqrt(v0.square()*v1.square()));
- const double d_alpha=alpha/(n+1);
const Point<spacedim> pm=0.5*(v0+v1);
const double h=std::sqrt(pm.square());
- double beta=0;
- unsigned int left_index=0, right_index=0;
- if ((n+1)%2==0)
- {
- // if the number of
- // parts is even insert
- // the midpoint
- left_index=(n-1)/2;
- right_index=left_index;
- points[left_index]=pm;
- --left_index;
- ++right_index;
- beta=d_alpha;
- }
- else
- {
- beta=(n/2+1)*d_alpha-alpha/2.;
- right_index=n/2;
- left_index=n/2-1;
- }
-
// n even: m=n/2,
// n odd: m=(n-1)/2
+ const std::vector<Point<1> > &line_points = this->get_line_support_points(n);
const unsigned int m=n/2;
- for (unsigned int i=0; i<m ; ++i, ++right_index, --left_index, beta+=d_alpha)
+ for (unsigned int i=0; i<m ; ++i)
{
- Assert(right_index<n, ExcInternalError());
- Assert(left_index+1>=1, ExcInternalError());
-
- double d=h*std::tan(beta);
- points[right_index]=pm+d/length*(v1-v0);
- points[left_index]=pm-d/length*(v1-v0);
+ const double beta = alpha * (line_points[i+1][0]-0.5);
+ const double d=h*std::tan(beta);
+ points[i]=pm+d/length*(v1-v0);
+ points[n-1-i]=pm-d/length*(v1-v0);
}
+ if ((n+1)%2==0)
+ // if the number of parts is even insert the midpoint
+ points[(n-1)/2]=pm;
- // project the points from the
- // straight line to the
- // HyperBallBoundary
+
+ // project the points from the straight line to the HyperBallBoundary
for (unsigned int i=0; i<n; ++i)
{
points[i] *= r / std::sqrt(points[i].square());
get_intermediate_points_on_line (const typename Triangulation<dim>::line_iterator &line,
std::vector<Point<dim> > &points) const
{
- // check whether center of object is
- // at x==0, since then it belongs
- // to the plane part of the
- // boundary
+ // check whether center of object is at x==0, since then it belongs to the
+ // plane part of the boundary
const Point<dim> line_center = line->center();
if (line_center(0) == this->center(0))
return StraightBoundary<dim>::get_intermediate_points_on_line (line, points);
points[0]=get_new_point_on_quad(quad);
else
{
- // check whether center of
- // object is at x==0, since
- // then it belongs to the plane
- // part of the boundary
+ // check whether center of object is at x==0, since then it belongs to
+ // the plane part of the boundary
const Point<dim> quad_center = quad->center();
if (quad_center(0) == this->center(0))
StraightBoundary<dim>::get_intermediate_points_on_quad (quad, points);
get_normals_at_vertices (const typename Triangulation<dim>::face_iterator &face,
typename Boundary<dim>::FaceVertexNormals &face_vertex_normals) const
{
- // check whether center of object is
- // at x==0, since then it belongs
- // to the plane part of the
- // boundary
+ // check whether center of object is at x==0, since then it belongs to the
+ // plane part of the boundary
const Point<dim> quad_center = face->center();
if (quad_center(0) == this->center(0))
StraightBoundary<dim>::get_normals_at_vertices (face, face_vertex_normals);
{
switch (dim)
{
- // in 2d, first check whether the two
- // end points of the line are on the
- // axis of symmetry. if so, then return
- // the mid point
+ // in 2d, first check whether the two end points of the line are on the
+ // axis of symmetry. if so, then return the mid point
case 2:
{
if ((line->vertex(0)(0) == this->center(0))
(line->vertex(1)(0) == this->center(0)))
return (line->vertex(0) + line->vertex(1))/2;
else
- // otherwise we are on the outer or
- // inner part of the shell. proceed
+ // otherwise we are on the outer or inner part of the shell. proceed
// as in the base class
return HyperShellBoundary<dim>::get_new_point_on_line (line);
}
- // in 3d, a line is a straight
- // line if it is on the symmetry
- // plane and if not both of its
- // end points are on either the
- // inner or outer sphere
+ // in 3d, a line is a straight line if it is on the symmetry plane and if
+ // not both of its end points are on either the inner or outer sphere
case 3:
{
- outer_radius) < 1e-12 * outer_radius))))
return (line->vertex(0) + line->vertex(1))/2;
else
- // otherwise we are on the outer or
- // inner part of the shell. proceed
+ // otherwise we are on the outer or inner part of the shell. proceed
// as in the base class
return HyperShellBoundary<dim>::get_new_point_on_line (line);
}
HalfHyperShellBoundary<dim>::
get_new_point_on_quad (const typename Triangulation<dim>::quad_iterator &quad) const
{
- // if this quad is on the symmetry plane,
- // take the center point and project it
- // outward to the same radius as the
- // centers of the two radial lines
+ // if this quad is on the symmetry plane, take the center point and project
+ // it outward to the same radius as the centers of the two radial lines
if ((quad->vertex(0)(0) == this->center(0)) &&
(quad->vertex(1)(0) == this->center(0)) &&
(quad->vertex(2)(0) == this->center(0)) &&
Assert (false, ExcInternalError());
}
- // otherwise we are on the outer or
- // inner part of the shell. proceed
- // as in the base class
+ // otherwise we are on the outer or inner part of the shell. proceed as in
+ // the base class
return HyperShellBoundary<dim>::get_new_point_on_quad (quad);
}
{
switch (dim)
{
- // in 2d, first check whether the two
- // end points of the line are on the
- // axis of symmetry. if so, then return
- // the mid point
+ // in 2d, first check whether the two end points of the line are on the
+ // axis of symmetry. if so, then return the mid point
case 2:
{
if ((line->vertex(0)(0) == this->center(0))
(line->vertex(1)(0) == this->center(0)))
StraightBoundary<dim>::get_intermediate_points_on_line (line, points);
else
- // otherwise we are on the outer or
- // inner part of the shell. proceed
+ // otherwise we are on the outer or inner part of the shell. proceed
// as in the base class
HyperShellBoundary<dim>::get_intermediate_points_on_line (line, points);
break;
}
- // in 3d, a line is a straight
- // line if it is on the symmetry
- // plane and if not both of its
- // end points are on either the
- // inner or outer sphere
+ // in 3d, a line is a straight line if it is on the symmetry plane and if
+ // not both of its end points are on either the inner or outer sphere
case 3:
{
if (((line->vertex(0)(0) == this->center(0))
- outer_radius) < 1e-12 * outer_radius))))
StraightBoundary<dim>::get_intermediate_points_on_line (line, points);
else
- // otherwise we are on the outer or
- // inner part of the shell. proceed
+ // otherwise we are on the outer or inner part of the shell. proceed
// as in the base class
HyperShellBoundary<dim>::get_intermediate_points_on_line (line, points);
{
Assert (dim < 3, ExcNotImplemented());
- // check whether center of object is
- // at x==0, since then it belongs
- // to the plane part of the
- // boundary
+ // check whether center of object is at x==0, since then it belongs to the
+ // plane part of the boundary
const Point<dim> quad_center = quad->center();
if (quad_center(0) == this->center(0))
StraightBoundary<dim>::get_intermediate_points_on_quad (quad, points);
Point<2> middle(0,0);
- //Take care for periodic conditions,
- //For instance phi0= 0, phi1= 3/2*Pi middle has to be 7/4*Pi not 3/4*Pi
- //This also works for -Pi/2 + Pi, middle is 5/4*Pi
+ //Take care for periodic conditions, For instance phi0= 0, phi1= 3/2*Pi
+ //middle has to be 7/4*Pi not 3/4*Pi. This also works for -Pi/2 + Pi, middle
+ //is 5/4*Pi
for (unsigned int i=0; i<2; i++)
if (std::abs(p0(i)-p1(i))> numbers::PI)
middle(i)=2*numbers::PI;
get_intermediate_points_on_line (const Triangulation<2, 3>::line_iterator &line,
std::vector< Point< 3 > > &points) const
{
- //Almost the same implementation as StraightBoundary<2,3>
+ //Almost the same implementation as StraightBoundary<2,3>
unsigned int npoints=points.size();
if (npoints==0) return;
offset[0]=0;
offset[1]=0;
- //Take care for periodic conditions & negative angles,
- //see get_new_point_on_line() above
- //Because we dont have a symmetric interpolation (just the middle) we need to
- //add 2*Pi to each almost zero and negative angles.
+ //Take care for periodic conditions & negative angles, see
+ //get_new_point_on_line() above. Because we dont have a symmetric
+ //interpolation (just the middle) we need to add 2*Pi to each almost zero
+ //and negative angles.
for (unsigned int i=0; i<2; i++)
for (unsigned int j=1; j<2; j++)
{
p[j](i)+=2*numbers::PI*offset[i];
- double dx=1.0/(npoints+1);
- double x=dx;
Point<2> target;
- for (unsigned int i=0; i<npoints; i++,x+=dx)
+ const std::vector<Point<1> > &line_points = this->get_line_support_points(npoints);
+ for (unsigned int i=0; i<npoints; i++)
{
+ const double x = line_points[i+1][0];
target= (1-x)*p[0] + x*p[1];
points[i]=get_real_coord(target);
}
// is n a square number
Assert(m*m==n, ExcInternalError());
- const double ds=1./(m+1);
- double y=ds;
-
Point<2> p[4];
for (unsigned int i=0; i<4; i++)
offset[0]=0;
offset[1]=0;
- //Take care for periodic conditions & negative angles,
- //see get_new_point_on_line() above
- //Because we dont have a symmetric interpolation (just the middle) we need to
- //add 2*Pi to each almost zero and negative angles.
+ //Take care for periodic conditions & negative angles, see
+ //get_new_point_on_line() above. Because we dont have a symmetric
+ //interpolation (just the middle) we need to add 2*Pi to each almost zero
+ //and negative angles.
for (unsigned int i=0; i<2; i++)
for (unsigned int j=1; j<4; j++)
{
if (p[j](i)<1.E-12 ) //Take care for periodic conditions & negative angles
p[j](i)+=2*numbers::PI*offset[i];
- for (unsigned int i=0; i<m; ++i, y+=ds)
+ const std::vector<Point<1> > &line_points = this->get_line_support_points(m);
+ for (unsigned int i=0; i<m; ++i)
{
- double x=ds;
- for (unsigned int j=0; j<m; ++j, x+=ds)
+ const double y=line_points[i+1][0];
+ for (unsigned int j=0; j<m; ++j)
{
+ const double x=line_points[j+1][0];
target=((1-x) * p[0] +
x * p[1]) * (1-y) +
((1-x) * p[2] +
DEAL::2d Jacobians:
-DEAL::0.9742 0.0664 -0.0314 0.5822
-DEAL::0.6722 0.0174 -0.0273 1.0693
+DEAL::0.9742 0.0664 -0.0313 0.5814
+DEAL::0.6710 0.0173 -0.0272 1.0693
DEAL::0.5858 0 0 0.5858
-DEAL::0.0314 -0.5822 0.9742 0.0664
-DEAL::-0.0273 1.0693 -0.6722 -0.0174
+DEAL::0.0313 -0.5814 0.9742 0.0664
+DEAL::-0.0272 1.0693 -0.6710 -0.0173
DEAL::
DEAL::2d inverse Jacobians:
-DEAL::1.0228 -0.1167 0.0552 1.7113
-DEAL::1.4866 -0.0242 0.0379 0.9345
+DEAL::1.0228 -0.1169 0.0550 1.7137
+DEAL::1.4893 -0.0241 0.0380 0.9345
DEAL::1.7071 0 0 1.7071
-DEAL::0.1167 1.0228 -1.7113 0.0552
-DEAL::-0.0242 -1.4866 0.9345 -0.0379
+DEAL::0.1169 1.0228 -1.7137 0.0550
+DEAL::-0.0241 -1.4893 0.9345 -0.0380
DEAL::
DEAL::3d Jacobians:
DEAL::0.4226 0 0 0 0.4226 0 0 0 0.4226
-DEAL::0.7166 0.0006 0.0633 0.0006 0.7163 -0.0221 -0.0143 0.0053 0.4812
-DEAL::0.0204 -0.5582 -0.0364 0.8017 0.0685 0.0034 0.0038 -0.1227 0.8075
-DEAL::0.8017 0.0685 0.0034 0.0038 -0.1227 0.8075 0.0204 -0.5582 -0.0364
-DEAL::0.6915 0.0121 0.0583 -0.0304 0.9077 -0.0021 -0.1531 -0.0021 0.9115
-DEAL::0.8017 0.0685 0.0034 -0.0204 0.5582 0.0364 0.0038 -0.1227 0.8075
-DEAL::-0.0304 0.9077 -0.0021 -0.6915 -0.0121 -0.0583 -0.1531 -0.0021 0.9115
+DEAL::0.7165 0.0006 0.0634 0.0006 0.7162 -0.0222 -0.0141 0.0052 0.4802
+DEAL::0.0201 -0.5556 -0.0358 0.8018 0.0685 0.0034 0.0038 -0.1227 0.8074
+DEAL::0.8018 0.0685 0.0034 0.0038 -0.1227 0.8074 0.0201 -0.5556 -0.0358
+DEAL::0.6859 0.0119 0.0572 -0.0304 0.9075 -0.0021 -0.1524 -0.0021 0.9114
+DEAL::0.8018 0.0685 0.0034 -0.0201 0.5556 0.0358 0.0038 -0.1227 0.8074
+DEAL::-0.0304 0.9075 -0.0021 -0.6859 -0.0119 -0.0572 -0.1524 -0.0021 0.9114
DEAL::
DEAL::3d inverse Jacobians:
DEAL::2.3660 0 0 0 2.3660 0 0 0 2.3660
-DEAL::1.3918 0.0002 -0.1831 0.0002 1.3957 0.0642 0.0414 -0.0155 2.0720
-DEAL::0.1523 1.2434 0.0016 -1.7683 0.0455 -0.0800 -0.2694 0.0011 1.2263
-DEAL::1.2434 0.0016 0.1523 0.0455 -0.0800 -1.7683 0.0011 1.2263 -0.2694
-DEAL::1.4251 -0.0193 -0.0912 0.0483 1.1010 -0.0005 0.2394 -0.0007 1.0818
-DEAL::1.2434 -0.1523 0.0016 0.0455 1.7683 -0.0800 0.0011 0.2694 1.2263
-DEAL::-0.0193 -1.4251 -0.0912 1.1010 -0.0483 -0.0005 -0.0007 -0.2394 1.0818
+DEAL::1.3921 0.0002 -0.1837 0.0001 1.3959 0.0646 0.0408 -0.0152 2.0764
+DEAL::0.1531 1.2433 0.0015 -1.7770 0.0449 -0.0790 -0.2707 0.0010 1.2265
+DEAL::1.2433 0.0015 0.1531 0.0449 -0.0790 -1.7770 0.0010 1.2265 -0.2707
+DEAL::1.4370 -0.0191 -0.0902 0.0487 1.1012 -0.0005 0.2404 -0.0006 1.0822
+DEAL::1.2433 -0.1531 0.0015 0.0449 1.7770 -0.0790 0.0010 0.2707 1.2265
+DEAL::-0.0191 -1.4370 -0.0902 1.1012 -0.0487 -0.0005 -0.0006 -0.2404 1.0822
DEAL::
DEAL:: area=6.8
DEAL::MappingFace2d-2-Q2
DEAL:MappingFace2d-2-Q2::1.0 1.1 -1.0 0.0 0.35
-DEAL:MappingFace2d-2-Q2::1.0 1.7 -1.0 -1.1e-16 0.65
+DEAL:MappingFace2d-2-Q2::1.0 1.7 -1.0 0.0 0.65
DEAL:MappingFace2d-2-Q2::1.0 2.3 -1.0 0.0 0.65
DEAL:MappingFace2d-2-Q2::1.0 2.9 -1.0 0.0 0.35
DEAL:MappingFace2d-2-Q2::
DEAL:MappingFace2d-2-Q2::
DEAL:MappingFace2d-2-Q2::
DEAL:MappingFace2d-2-Q2::1.3 1.0 0.0 -1.0 0.70
-DEAL:MappingFace2d-2-Q2::2.3 1.0 -5.6e-17 -1.0 1.3
+DEAL:MappingFace2d-2-Q2::2.3 1.0 0.0 -1.0 1.3
DEAL:MappingFace2d-2-Q2::3.7 1.0 0.0 -1.0 1.3
DEAL:MappingFace2d-2-Q2::4.7 1.0 0.0 -1.0 0.70
DEAL:MappingFace2d-2-Q2::
DEAL:MappingFace2d-2-Q2::
DEAL:MappingFace2d-2-Q2::
DEAL::MappingSubface2d-2-Q2
-DEAL:MappingSubface2d-2-Q2::1.0 1.1 -1.0 -2.2e-16
+DEAL:MappingSubface2d-2-Q2::1.0 1.1 -1.0 0.0
DEAL:MappingSubface2d-2-Q2::1.0 1.3 -1.0 0.0
DEAL:MappingSubface2d-2-Q2::1.0 1.7 -1.0 0.0
-DEAL:MappingSubface2d-2-Q2::1.0 1.9 -1.0 -1.1e-16
+DEAL:MappingSubface2d-2-Q2::1.0 1.9 -1.0 0.0
DEAL:MappingSubface2d-2-Q2::
DEAL:MappingSubface2d-2-Q2::
DEAL:MappingSubface2d-2-Q2::1.0 2.1 -1.0 0.0
DEAL:MappingSubface2d-2-Q2::3.1 3.0 0.13 0.99
DEAL:MappingSubface2d-2-Q2::
DEAL:MappingSubface2d-2-Q2::
-DEAL:MappingSubface2d-2-Q2::1.1 1.0 -1.1e-16 -1.0
+DEAL:MappingSubface2d-2-Q2::1.1 1.0 0.0 -1.0
DEAL:MappingSubface2d-2-Q2::1.7 1.0 0.0 -1.0
DEAL:MappingSubface2d-2-Q2::2.3 1.0 0.0 -1.0
-DEAL:MappingSubface2d-2-Q2::2.9 1.0 -5.6e-17 -1.0
+DEAL:MappingSubface2d-2-Q2::2.9 1.0 0.0 -1.0
DEAL:MappingSubface2d-2-Q2::
DEAL:MappingSubface2d-2-Q2::
DEAL:MappingSubface2d-2-Q2::3.1 1.0 0.0 -1.0
DEAL:Mapping2d-2-Q3::2.0 1.6 7.1
DEAL:Mapping2d-2-Q3::2.5 1.6 7.4
DEAL:Mapping2d-2-Q3::3.0 1.6 8.1
-DEAL:Mapping2d-2-Q3::3.6 1.7 9.4
+DEAL:Mapping2d-2-Q3::3.6 1.7 9.3
DEAL:Mapping2d-2-Q3::4.2 1.8 11.
DEAL:Mapping2d-2-Q3::4.8 1.9 14.
DEAL:Mapping2d-2-Q3::
-DEAL:Mapping2d-2-Q3::1.0 1.9 6.7
+DEAL:Mapping2d-2-Q3::1.0 1.9 6.6
DEAL:Mapping2d-2-Q3::1.5 1.9 6.4
DEAL:Mapping2d-2-Q3::1.9 1.9 6.3
DEAL:Mapping2d-2-Q3::2.4 1.9 6.6
DEAL:Mapping2d-2-Q3::2.3 2.1 5.6
DEAL:Mapping2d-2-Q3::2.7 2.2 6.1
DEAL:Mapping2d-2-Q3::3.2 2.3 7.1
-DEAL:Mapping2d-2-Q3::3.7 2.4 8.8
+DEAL:Mapping2d-2-Q3::3.7 2.4 8.7
DEAL:Mapping2d-2-Q3::4.2 2.6 11.
DEAL:Mapping2d-2-Q3::
DEAL:Mapping2d-2-Q3::1.0 2.4 5.3
DEAL:Mapping2d-2-Q3::1.3 2.7 3.8
DEAL:Mapping2d-2-Q3::1.7 2.6 3.4
DEAL:Mapping2d-2-Q3::2.0 2.6 3.3
-DEAL:Mapping2d-2-Q3::2.4 2.6 3.4
+DEAL:Mapping2d-2-Q3::2.4 2.6 3.5
DEAL:Mapping2d-2-Q3::2.7 2.7 3.9
DEAL:Mapping2d-2-Q3::3.1 2.8 4.7
DEAL:Mapping2d-2-Q3::3.4 3.0 6.1
DEAL:Mapping2d-2-Q3::
DEAL:: area=6.8
DEAL::MappingFace2d-2-Q3
-DEAL:MappingFace2d-2-Q3::1.0 1.1 -1.0 4.4e-16 0.35
+DEAL:MappingFace2d-2-Q3::1.0 1.1 -1.0 -7.8e-16 0.35
DEAL:MappingFace2d-2-Q3::1.0 1.7 -1.0 0.0 0.65
-DEAL:MappingFace2d-2-Q3::1.0 2.3 -1.0 -6.7e-16 0.65
+DEAL:MappingFace2d-2-Q3::1.0 2.3 -1.0 5.6e-17 0.65
DEAL:MappingFace2d-2-Q3::1.0 2.9 -1.0 0.0 0.35
DEAL:MappingFace2d-2-Q3::
DEAL:MappingFace2d-2-Q3::
-DEAL:MappingFace2d-2-Q3::5.0 1.2 1.0 0.098 0.55
+DEAL:MappingFace2d-2-Q3::5.0 1.2 0.99 0.10 0.55
DEAL:MappingFace2d-2-Q3::4.7 2.0 0.86 0.51 1.0
DEAL:MappingFace2d-2-Q3::4.0 2.7 0.51 0.86 1.0
-DEAL:MappingFace2d-2-Q3::3.2 3.0 0.098 1.0 0.55
+DEAL:MappingFace2d-2-Q3::3.2 3.0 0.10 0.99 0.55
DEAL:MappingFace2d-2-Q3::
DEAL:MappingFace2d-2-Q3::
-DEAL:MappingFace2d-2-Q3::1.3 1.0 2.2e-16 -1.0 0.70
+DEAL:MappingFace2d-2-Q3::1.3 1.0 -3.9e-16 -1.0 0.70
DEAL:MappingFace2d-2-Q3::2.3 1.0 0.0 -1.0 1.3
-DEAL:MappingFace2d-2-Q3::3.7 1.0 -3.3e-16 -1.0 1.3
+DEAL:MappingFace2d-2-Q3::3.7 1.0 2.8e-17 -1.0 1.3
DEAL:MappingFace2d-2-Q3::4.7 1.0 0.0 -1.0 0.70
DEAL:MappingFace2d-2-Q3::
DEAL:MappingFace2d-2-Q3::
DEAL:MappingFace2d-2-Q3::
DEAL:MappingFace2d-2-Q3::
DEAL::MappingSubface2d-2-Q3
-DEAL:MappingSubface2d-2-Q3::1.0 1.1 -1.0 -4.4e-16
-DEAL:MappingSubface2d-2-Q3::1.0 1.3 -1.0 -2.2e-16
-DEAL:MappingSubface2d-2-Q3::1.0 1.7 -1.0 -4.4e-16
-DEAL:MappingSubface2d-2-Q3::1.0 1.9 -1.0 2.2e-16
+DEAL:MappingSubface2d-2-Q3::1.0 1.1 -1.0 -8.9e-16
+DEAL:MappingSubface2d-2-Q3::1.0 1.3 -1.0 -5.0e-16
+DEAL:MappingSubface2d-2-Q3::1.0 1.7 -1.0 0.0
+DEAL:MappingSubface2d-2-Q3::1.0 1.9 -1.0 -2.2e-16
DEAL:MappingSubface2d-2-Q3::
DEAL:MappingSubface2d-2-Q3::
-DEAL:MappingSubface2d-2-Q3::1.0 2.1 -1.0 -8.9e-16
-DEAL:MappingSubface2d-2-Q3::1.0 2.3 -1.0 -1.1e-16
+DEAL:MappingSubface2d-2-Q3::1.0 2.1 -1.0 2.2e-16
+DEAL:MappingSubface2d-2-Q3::1.0 2.3 -1.0 0.0
DEAL:MappingSubface2d-2-Q3::1.0 2.7 -1.0 0.0
-DEAL:MappingSubface2d-2-Q3::1.0 2.9 -1.0 -4.4e-16
+DEAL:MappingSubface2d-2-Q3::1.0 2.9 -1.0 4.4e-16
DEAL:MappingSubface2d-2-Q3::
DEAL:MappingSubface2d-2-Q3::
-DEAL:MappingSubface2d-2-Q3::5.0 1.1 1.0 0.036
-DEAL:MappingSubface2d-2-Q3::4.9 1.5 0.97 0.26
+DEAL:MappingSubface2d-2-Q3::5.0 1.1 1.0 0.039
+DEAL:MappingSubface2d-2-Q3::4.9 1.5 0.96 0.26
DEAL:MappingSubface2d-2-Q3::4.7 2.0 0.86 0.51
DEAL:MappingSubface2d-2-Q3::4.5 2.3 0.74 0.67
DEAL:MappingSubface2d-2-Q3::
DEAL:MappingSubface2d-2-Q3::
DEAL:MappingSubface2d-2-Q3::4.3 2.5 0.67 0.74
DEAL:MappingSubface2d-2-Q3::4.0 2.7 0.51 0.86
-DEAL:MappingSubface2d-2-Q3::3.5 2.9 0.26 0.97
-DEAL:MappingSubface2d-2-Q3::3.1 3.0 0.036 1.0
+DEAL:MappingSubface2d-2-Q3::3.5 2.9 0.26 0.96
+DEAL:MappingSubface2d-2-Q3::3.1 3.0 0.039 1.0
DEAL:MappingSubface2d-2-Q3::
DEAL:MappingSubface2d-2-Q3::
-DEAL:MappingSubface2d-2-Q3::1.1 1.0 -2.2e-16 -1.0
-DEAL:MappingSubface2d-2-Q3::1.7 1.0 -1.1e-16 -1.0
-DEAL:MappingSubface2d-2-Q3::2.3 1.0 -2.2e-16 -1.0
-DEAL:MappingSubface2d-2-Q3::2.9 1.0 1.1e-16 -1.0
+DEAL:MappingSubface2d-2-Q3::1.1 1.0 -4.4e-16 -1.0
+DEAL:MappingSubface2d-2-Q3::1.7 1.0 -2.5e-16 -1.0
+DEAL:MappingSubface2d-2-Q3::2.3 1.0 0.0 -1.0
+DEAL:MappingSubface2d-2-Q3::2.9 1.0 -1.1e-16 -1.0
DEAL:MappingSubface2d-2-Q3::
DEAL:MappingSubface2d-2-Q3::
-DEAL:MappingSubface2d-2-Q3::3.1 1.0 -4.4e-16 -1.0
-DEAL:MappingSubface2d-2-Q3::3.7 1.0 -5.6e-17 -1.0
+DEAL:MappingSubface2d-2-Q3::3.1 1.0 1.1e-16 -1.0
+DEAL:MappingSubface2d-2-Q3::3.7 1.0 0.0 -1.0
DEAL:MappingSubface2d-2-Q3::4.3 1.0 0.0 -1.0
-DEAL:MappingSubface2d-2-Q3::4.9 1.0 -2.2e-16 -1.0
+DEAL:MappingSubface2d-2-Q3::4.9 1.0 2.2e-16 -1.0
DEAL:MappingSubface2d-2-Q3::
DEAL:MappingSubface2d-2-Q3::
DEAL:MappingSubface2d-2-Q3::1.1 3.0 0.42 0.91
4 162 3.14159639 1.1900e-06 5.92
16 578 3.14159271 1.8860e-08 5.98
64 2178 3.14159265 2.9574e-10 5.99
-256 8450 3.14159265 4.6257e-12 6.00
+256 8450 3.14159265 4.6356e-12 6.00
DEAL::
// initial guess for the inverse Q4 mapping by running an inverse Q1
// mapping, and we failed to find a solution in the latter.
//
+// todo: this test used to fail before r30450 but the new mapping support
+// points seem to fix this
+//
// a redux of this test is in the mapping_real_to_unit_q4_sphere_*
// tests
projects / dealii.git / commitdiff