* <h4>Behavior along curved boundaries and with different manifolds</h4>
*
* As described above, one often only knows a manifold description of a
- * surface but not the interior of the computational domain. Thus, a
+ * surface but not the interior of the computational domain. In such a case, a
* StraightBoundary object will be assigned to the interior entities that
- * describes a usual planar coordinate system where the additional points of
- * this class are placed exactly according to a bi-/trilinear mapping. When
- * combined with a non-flat manifold on the boundary, for example a circle,
- * the two manifold descriptions are in general incompatible. For example, a
- * StraightBoundary would put an interior point located at some small distance
+ * describes a usual planar coordinate system where the additional points for
+ * the higher order mapping are placed exactly according to a bi-/trilinear
+ * mapping. When combined with a non-flat manifold on the boundary, for
+ * example a circle, the two manifold descriptions are in general
+ * incompatible. For example, a StraightBoundary defined solely through the
+ * cell's vertices would put an interior point located at some small distance
* epsilon away from the boundary along a flat line and thus in general
* outside the concave part of a circle. If the polynomial degree of
* MappingQGeneric is sufficiently high, the transformation from the reference
* boundary.
*
* In order to avoid this situation, this class applies a smoothing on cells
- * adjacent to the boundary by using so-called Laplace smoothing. In the
- * algorithm computing new points, all the entities of the cells are passed
- * through hierarchically, starting from the lines to the quads and finally
- * hexes. The elements higher up in the hierarchy that sit on a
- * StraightBoundary will then get their points interpolated from all the
- * surrounding points and not just the corner points. If only a line is
- * associated to a curved boundary but the adjacent quad is on a flat
- * manifold, the points inside the quad will be computed according to the
- * deformed line and thus always result in a well-defined transformation. If
- * this smoothing is undesired, the optional argument @p
- * smooth_support_points_on_flat_manifold can be used to disable the
- * smoothing, placing the additional points strictly according to the
- * manifold.
+ * adjacent to the boundary by using so-called Laplace smoothing by
+ * default. In the algorithm that computing additional points, the
+ * compute_mapping_support_points() method, all the entities of the cells are
+ * passed through hierarchically, starting from the lines to the quads and
+ * finally hexes. The elements higher up in the hierarchy that sit on the
+ * boundary will then get their points interpolated from all the surrounding
+ * points and not just the corner points. If only a line is assigned a curved
+ * boundary but the adjacent quad is on a flat manifold, the points inside the
+ * quad will be computed according to the deformed line and thus always result
+ * in a well-defined transformation. This smoothing can be disabled by setting
+ * the optional argument @p smooth_support_points to false, placing the
+ * additional points strictly according to the manifold. This is usually the
+ * most efficient choice in case different manifolds are present that are
+ * compatible with each other.
*
- * While this smoothing approach works well for low and medium convergence
- * orders up to approximately three to four, this mechanism has an inherent
- * shortcoming because it switches from a curved manifold to a flat manifold
- * within one layer of elements. This will cause the Jacobian transformation
- * to have jumps between the first and second element layer that can reduce
- * the order of convergence. For example, the convergence rates for solving
- * the Laplacian on a circle where only the boundary is deformed and the above
- * mesh smoothing algorithm is applied will typically not exceed 3.5, even for
+ * While the smoothing approach works well for filling holes or avoiding
+ * inversions with low and medium convergence orders up to approximately three
+ * to four, there is nonetheless an inherent shortcoming because of a
+ * discontinuous mapping that switches from a curved manifold to a flat
+ * manifold within one layer of elements. This will cause the Jacobian
+ * transformation to have jumps between the first and second element layer
+ * that can reduce the order of convergence. For example, the convergence
+ * rates for solving the Laplacian on a circle where only the boundary is
+ * deformed and the above mesh smoothing algorithm is applied will typically
+ * not exceed 3.5 (or 3 in the elements adjacent to the boundary), even for
* fourth or fifth degree polynomials. In such a case, the curved manifold
* needs to be switched to a flat manifold in a smooth way that does not
* depend on the mesh size and eventuell covers a whole layer of cells.
* equal to @p spacedim, i.e., if we are considering meshes on surfaces
* embedded into higher dimensional spaces.
*
- * The optional parameter @p smooth_support_points_on_flat_manifold controls
- * whether smoothing adjancent to StraightBoundary objects according to the
+ * The optional parameter @p smooth_support_points controls whether
+ * smoothing on objects where different manifolds meet according to the
* general class description should be enabled (default) or not.
*/
MappingQ (const unsigned int polynomial_degree,
const bool use_mapping_q_on_all_cells = false,
- const bool smooth_support_points_on_flat_manifold = true);
+ const bool smooth_support_points = true);
/**
* Copy constructor.
* <h4>Behavior along curved boundaries and with different manifolds</h4>
*
* As described above, one often only knows a manifold description of a
- * surface but not the interior of the computational domain. Thus, a
+ * surface but not the interior of the computational domain. In such a case, a
* StraightBoundary object will be assigned to the interior entities that
- * describes a usual planar coordinate system where the additional points of
- * this class are placed exactly according to a bi-/trilinear mapping. When
- * combined with a non-flat manifold on the boundary, for example a circle,
- * the two manifold descriptions are in general incompatible. For example, a
- * StraightBoundary would put an interior point located at some small distance
+ * describes a usual planar coordinate system where the additional points for
+ * the higher order mapping are placed exactly according to a bi-/trilinear
+ * mapping. When combined with a non-flat manifold on the boundary, for
+ * example a circle, the two manifold descriptions are in general
+ * incompatible. For example, a StraightBoundary defined solely through the
+ * cell's vertices would put an interior point located at some small distance
* epsilon away from the boundary along a flat line and thus in general
* outside the concave part of a circle. If the polynomial degree of
* MappingQGeneric is sufficiently high, the transformation from the reference
* boundary.
*
* In order to avoid this situation, this class applies a smoothing on cells
- * adjacent to the boundary by using so-called Laplace smoothing. In the
- * algorithm computing new points, all the entities of the cells are passed
- * through hierarchically, starting from the lines to the quads and finally
- * hexes. The elements higher up in the hierarchy that sit on a
- * StraightBoundary will then get their points interpolated from all the
- * surrounding points and not just the corner points. If only a line is
- * associated to a curved boundary but the adjacent quad is on a flat
- * manifold, the points inside the quad will be computed according to the
- * deformed line and thus always result in a well-defined transformation. If
- * this smoothing is undesired, the optional argument @p
- * smooth_support_points_on_flat_manifold can be used to disable the
- * smoothing, placing the additional points strictly according to the
- * manifold.
+ * adjacent to the boundary by using so-called Laplace smoothing by
+ * default. In the algorithm that computing additional points, the
+ * compute_mapping_support_points() method, all the entities of the cells are
+ * passed through hierarchically, starting from the lines to the quads and
+ * finally hexes. The elements higher up in the hierarchy that sit on the
+ * boundary will then get their points interpolated from all the surrounding
+ * points and not just the corner points. If only a line is assigned a curved
+ * boundary but the adjacent quad is on a flat manifold, the points inside the
+ * quad will be computed according to the deformed line and thus always result
+ * in a well-defined transformation. This smoothing can be disabled by setting
+ * the optional argument @p smooth_support_points to false, placing the
+ * additional points strictly according to the manifold. This is usually the
+ * most efficient choice in case different manifolds are present that are
+ * compatible with each other.
*
- * While this smoothing approach works well for low and medium convergence
- * orders up to approximately three to four, this mechanism has an inherent
- * shortcoming because it switches from a curved manifold to a flat manifold
- * within one layer of elements. This will cause the Jacobian transformation
- * to have jumps between the first and second element layer that can reduce
- * the order of convergence. For example, the convergence rates for solving
- * the Laplacian on a circle where only the boundary is deformed and the above
- * mesh smoothing algorithm is applied will typically not exceed 3.5, even for
+ * While the smoothing approach works well for filling holes or avoiding
+ * inversions with low and medium convergence orders up to approximately three
+ * to four, there is nonetheless an inherent shortcoming because of a
+ * discontinuous mapping that switches from a curved manifold to a flat
+ * manifold within one layer of elements. This will cause the Jacobian
+ * transformation to have jumps between the first and second element layer
+ * that can reduce the order of convergence. For example, the convergence
+ * rates for solving the Laplacian on a circle where only the boundary is
+ * deformed and the above mesh smoothing algorithm is applied will typically
+ * not exceed 3.5 (or 3 in the elements adjacent to the boundary), even for
* fourth or fifth degree polynomials. In such a case, the curved manifold
* needs to be switched to a flat manifold in a smooth way that does not
* depend on the mesh size and eventuell covers a whole layer of cells.
*
- * @author Wolfgang Bangerth, 2015
+ * @author Wolfgang Bangerth, 2015, Martin Kronbichler, 2017
*/
template <int dim, int spacedim=dim>
class MappingQGeneric : public Mapping<dim,spacedim>
* polynomials that are used to map cells from the reference to the real
* cell.
*
- * The optional parameter @p smooth_support_points_on_flat_manifold controls
- * whether smoothing adjancent to StraightBoundary objects according to the
- * general class description should be enabled (default) or not.
+ * The optional parameter @p smooth_support_points controls whether
+ * smoothing according to the general class description should be enabled
+ * (default) or not.
*/
MappingQGeneric (const unsigned int polynomial_degree,
- const bool smooth_support_points_on_flat_manifold = true);
+ const bool smooth_support_points = true);
/**
* Copy constructor.
*/
const unsigned int n_shape_functions;
- /*
- * The default line support points. Is used in when the shape function
- * values are computed.
- *
- * The number of quadrature points depends on the degree of this
- * class, and it matches the number of degrees of freedom of an
- * FE_Q<1>(this->degree).
- */
- QGaussLobatto<1> line_support_points;
-
/**
* Tensors of covariant transformation at each of the quadrature points.
* The matrix stored is the Jacobian * G^{-1}, where G = Jacobian^{t} *
* Stores whether we want to smooth the placement of interior points on flat
* manifolds.
*/
- const bool smooth_support_points_on_flat_manifold;
+ const bool smooth_support_points;
/*
* The default line support points. These are used when computing
const std_cxx11::unique_ptr<FE_Q<dim> > fe_q;
/**
- * A table of weights by which we multiply the locations of the support
- * points on the perimeter of a quad to get the location of interior support
- * points.
+ * A vector of tables of weights by which we multiply the locations of the
+ * support points on the perimeter of an object (line, quad, hex) to get the
+ * location of interior support points.
*
- * Sizes: support_point_weights_on_quad.size()= number of inner
- * unit_support_points support_point_weights_on_quad[i].size()= number of
- * outer unit_support_points, i.e. unit_support_points on the boundary of
- * the quad
+ * Access into this table is by @p [structdim-1], i.e., use 0 to access the
+ * support point weights on a line (i.e., the interior points of the
+ * GaussLobatto quadrature), use 1 to access the support point weights from
+ * to perimeter to the interior of a quad, and use 2 to access the support
+ * point weights from the perimeter to the interior of a hex.
*
- * For the definition of this vector see equation (8) of the `mapping'
+ * The table itself contains as many columns as there are surrounding points
+ * to a particular object (2 for a vertex, <code>4 + 4*(degree-1)</code> for
+ * a quad, <code>8 + 12*(degree-1) + 6*(degree-1)*(degree-1)</code> for a
+ * hex) and as many rows as there are strictly interior points.
+ *
+ * For the definition of this table see equation (8) of the `mapping'
* report.
*/
- Table<2,double> support_point_weights_on_quad;
+ std::vector<Table<2,double> > support_point_weights_perimeter_to_interior;
/**
- * A table of weights by which we multiply the locations of the support
- * points on the perimeter of a hex to get the location of interior support
- * points.
+ * A vector of tables of weights by which we multiply the locations of the
+ * vertex points of an object (line, quad, hex) to get the location of
+ * interior support points.
*
- * For the definition of this vector see equation (8) of the `mapping'
- * report.
+ * As opposed to @p support_point_weights_perimeter_to_interior, this table
+ * takes only the vertex points into account and not intermediate points
+ * inside the objects. Thus, the content of this table simply corresponds to
+ * the evaluation of the linear shape functions in the dimension of the
+ * object.
+ *
+ * Access into the vector of tables is by @p [structdim-1], i.e., use 0 to
+ * access the support point weights on a line (i.e., the interior points of
+ * the GaussLobatto quadrature), use 1 to access the support point weights
+ * from to perimeter to the interior of a quad, and use 2 to access the
+ * support point weights from the perimeter to the interior of a hex.
+ *
+ * The table itself contains as many columns as there are vertices on a
+ * particular object (2 for a vertex, 4 for a quad, 8 for a hex) and as many
+ * rows as there are strictly interior points.
+ */
+ std::vector<Table<2,double> > support_point_weights_interior;
+
+ /**
+ * A table of weights by which we multiply the locations of the vertex
+ * points of the cell to get the location of all additional support points,
+ * both on lines, quads, and hexes (as appropriate). This data structure is
+ * used when we fill all support points at once, which is the case if the
+ * same manifold is attached to all sub-entities of a cell. This way, we can
+ * avoid some of the overhead in transforming data for mappings.
+ *
+ * The table has as many rows as there are vertices to the cell (2 in 1D, 4
+ * in 2D, 8 in 3D), and as many rows as there are additional support points
+ * in the mapping, i.e., <code>(degree+1)^dim - 2^dim</code>.
*/
- Table<2,double> support_point_weights_on_hex;
+ Table<2,double> support_point_weights_cell;
/**
* Return the locations of support points for the mapping. For example, for
// Construct the tensor product polynomials used as shape functions for the
// Qp mapping of cells at the boundary.
+ QGaussLobatto<1> line_support_points(data.polynomial_degree+1);
const TensorProductPolynomials<dim>
- tensor_pols (Polynomials::generate_complete_Lagrange_basis(data.line_support_points.get_points()));
+ tensor_pols (Polynomials::generate_complete_Lagrange_basis(line_support_points.get_points()));
Assert (n_shape_functions==tensor_pols.n(),
ExcInternalError());
MappingQGeneric<dim,spacedim>::InternalData::InternalData (const unsigned int polynomial_degree)
:
polynomial_degree (polynomial_degree),
- n_shape_functions (Utilities::fixed_power<dim>(polynomial_degree+1)),
- line_support_points (polynomial_degree + 1)
+ n_shape_functions (Utilities::fixed_power<dim>(polynomial_degree+1))
{}
* For the definition of the @p support_point_weights_on_quad please refer to
* equation (8) of the `mapping' report.
*/
- template<int dim>
Table<2,double>
compute_support_point_weights_on_quad(const unsigned int polynomial_degree)
{
Table<2,double> loqvs;
- // in 1d, there are no quads, so return an empty object
- if (dim == 1)
- return loqvs;
-
// we are asked to compute weights for interior support points, but
// there are no interior points if degree==1
if (polynomial_degree == 1)
* For the definition of the @p support_point_weights_on_hex please refer to
* equation (8) of the `mapping' report.
*/
- template <int dim>
Table<2,double>
compute_support_point_weights_on_hex(const unsigned int polynomial_degree)
{
Table<2,double> lohvs;
- // in 1d and 2d, there are no hexes, so return an empty object
- if (dim < 3)
- return lohvs;
-
// we are asked to compute weights for interior support points, but
// there are no interior points if degree==1
if (polynomial_degree == 1)
return lohvs;
- const unsigned int n_inner = Utilities::fixed_power<dim>(polynomial_degree-1);
+ const unsigned int n_inner = Utilities::fixed_power<3>(polynomial_degree-1);
const unsigned int n_outer = 8+12*(polynomial_degree-1)+6*(polynomial_degree-1)*(polynomial_degree-1);
// first check whether we have precomputed the values for some polynomial
else
{
// not precomputed, then do so now
- lohvs = compute_laplace_vector<dim>(polynomial_degree);
+ lohvs = compute_laplace_vector<3>(polynomial_degree);
}
// the sum of weights of the points at the outer rim should be one. check
return lohvs;
}
+
+ /**
+ * This function collects the output of
+ * compute_support_point_weights_on_{quad,hex} in a single data structure.
+ */
+ std::vector<Table<2,double> >
+ compute_support_point_weights_perimeter_to_interior(const unsigned int polynomial_degree,
+ const unsigned int dim)
+ {
+ Assert(dim > 0 && dim <= 3, ExcImpossibleInDim(dim));
+ std::vector<Table<2,double> > output(dim);
+ if (polynomial_degree <= 1)
+ return output;
+
+ // fill the 1D interior weights
+ QGaussLobatto<1> quadrature(polynomial_degree+1);
+ output[0].reinit(polynomial_degree-1, GeometryInfo<1>::vertices_per_cell);
+ for (unsigned int q=0; q<polynomial_degree-1; ++q)
+ for (unsigned int i=0; i<GeometryInfo<1>::vertices_per_cell; ++i)
+ output[0](q,i) = GeometryInfo<1>::d_linear_shape_function(quadrature.point(q+1),
+ i);
+
+ if (dim > 1)
+ output[1] = compute_support_point_weights_on_quad(polynomial_degree);
+
+ if (dim > 2)
+ output[2] = compute_support_point_weights_on_hex(polynomial_degree);
+
+ return output;
+ }
+
+ /**
+ * Collects all interior points for the various dimensions.
+ */
+ std::vector<Table<2,double> >
+ compute_support_point_weights_interior(const unsigned int polynomial_degree,
+ const unsigned int dim)
+ {
+ Assert(dim > 0 && dim <= 3, ExcImpossibleInDim(dim));
+ std::vector<Table<2,double> > output(dim);
+ if (polynomial_degree <= 1)
+ return output;
+
+ // fill the 1D interior weights
+ QGaussLobatto<1> quadrature(polynomial_degree+1);
+ output[0].reinit(polynomial_degree-1, GeometryInfo<1>::vertices_per_cell);
+ for (unsigned int q=0; q<polynomial_degree-1; ++q)
+ for (unsigned int i=0; i<GeometryInfo<1>::vertices_per_cell; ++i)
+ output[0](q,i) = GeometryInfo<1>::d_linear_shape_function(quadrature.point(q+1),
+ i);
+
+ // fill the 2D interior points
+ if (dim > 1)
+ {
+ output[1].reinit(Utilities::fixed_power<2>(polynomial_degree-1),
+ GeometryInfo<2>::vertices_per_cell);
+ for (unsigned int q=0, q2=0; q2<polynomial_degree-1; ++q2)
+ for (unsigned int q1=0; q1<polynomial_degree-1; ++q1, ++q)
+ {
+ Point<2> point(quadrature.point(q1+1)[0],
+ quadrature.point(q2+1)[0]);
+ for (unsigned int i=0; i<GeometryInfo<2>::vertices_per_cell; ++i)
+ output[1](q,i) = GeometryInfo<2>::d_linear_shape_function(point, i);
+ }
+ }
+
+ // fill the 3D interior points
+ if (dim > 2)
+ {
+ output[2].reinit(Utilities::fixed_power<3>(polynomial_degree-1),
+ GeometryInfo<3>::vertices_per_cell);
+ for (unsigned int q=0, q3=0; q3<polynomial_degree-1; ++q3)
+ for (unsigned int q2=0; q2<polynomial_degree-1; ++q2)
+ for (unsigned int q1=0; q1<polynomial_degree-1; ++q1, ++q)
+ {
+ Point<3> point(quadrature.point(q1+1)[0],
+ quadrature.point(q2+1)[0],
+ quadrature.point(q3+1)[0]);
+ for (unsigned int i=0; i<GeometryInfo<3>::vertices_per_cell; ++i)
+ output[2](q,i) = GeometryInfo<3>::d_linear_shape_function(point, i);
+ }
+ }
+
+ return output;
+ }
+
+ /**
+ * Collects all interior points for the various dimensions.
+ */
+ template <int dim>
+ Table<2,double>
+ compute_support_point_weights_cell(const unsigned int polynomial_degree)
+ {
+ Assert(dim > 0 && dim <= 3, ExcImpossibleInDim(dim));
+ if (polynomial_degree <= 1)
+ return Table<2,double>();
+
+ FE_Q<dim> fe(polynomial_degree);
+ std::vector<Point<dim> > unit_support_points = fe.get_unit_support_points();
+
+ Table<2,double> output(fe.dofs_per_cell - GeometryInfo<dim>::vertices_per_cell,
+ GeometryInfo<dim>::vertices_per_cell);
+ for (unsigned int q=0; q<output.size(0); ++q)
+ for (unsigned int i=0; i<GeometryInfo<dim>::vertices_per_cell; ++i)
+ output(q,i) = GeometryInfo<dim>::d_linear_shape_function(unit_support_points[q+GeometryInfo<dim>::vertices_per_cell],
+ i);
+
+ return output;
+ }
}
template<int dim, int spacedim>
MappingQGeneric<dim,spacedim>::MappingQGeneric (const unsigned int p,
- const bool smooth_support_points_on_flat_manifold)
+ const bool smooth_support_points)
:
polynomial_degree(p),
- smooth_support_points_on_flat_manifold(smooth_support_points_on_flat_manifold),
+ smooth_support_points(smooth_support_points),
line_support_points(this->polynomial_degree+1),
fe_q(dim == 3 ? new FE_Q<dim>(this->polynomial_degree) : 0),
- support_point_weights_on_quad (compute_support_point_weights_on_quad<dim>(this->polynomial_degree)),
- support_point_weights_on_hex (compute_support_point_weights_on_hex<dim>(this->polynomial_degree))
+ support_point_weights_perimeter_to_interior (compute_support_point_weights_perimeter_to_interior(this->polynomial_degree, dim)),
+ support_point_weights_interior (compute_support_point_weights_interior(this->polynomial_degree, dim)),
+ support_point_weights_cell (compute_support_point_weights_cell<dim>(this->polynomial_degree))
{
Assert (p >= 1, ExcMessage ("It only makes sense to create polynomial mappings "
"with a polynomial degree greater or equal to one."));
MappingQGeneric<dim,spacedim>::MappingQGeneric (const MappingQGeneric<dim,spacedim> &mapping)
:
polynomial_degree(mapping.polynomial_degree),
- smooth_support_points_on_flat_manifold(mapping.smooth_support_points_on_flat_manifold),
+ smooth_support_points(mapping.smooth_support_points),
line_support_points(mapping.line_support_points),
fe_q(dim == 3 ? new FE_Q<dim>(*mapping.fe_q) : 0),
- support_point_weights_on_quad (mapping.support_point_weights_on_quad),
- support_point_weights_on_hex (mapping.support_point_weights_on_hex)
+ support_point_weights_perimeter_to_interior (mapping.support_point_weights_perimeter_to_interior),
+ support_point_weights_interior (mapping.support_point_weights_interior),
+ support_point_weights_cell (mapping.support_point_weights_cell)
{}
namespace
{
- /**
- * Ask the manifold descriptor to return intermediate points on lines or
- * faces. The function needs to return one or multiple points (depending on
- * the number of elements in the output vector @p points that lie inside a
- * line, quad or hex). Whether it is a line, quad or hex doesn't really
- * matter to this function but it can be inferred from the number of input
- * points in the @p surrounding_points vector.
- */
- template<int dim, int spacedim>
- void
- get_intermediate_points (const Manifold<dim, spacedim> &manifold,
- const QGaussLobatto<1> &line_support_points,
- const std::vector<Point<spacedim> > &surrounding_points,
- std::vector<Point<spacedim> > &points)
- {
- Assert(surrounding_points.size() >= 2, ExcMessage("At least 2 surrounding points are required"));
- const unsigned int n=points.size();
- Assert(n>0, ExcMessage("You can't ask for 0 intermediate points."));
- std::vector<double> w(surrounding_points.size());
-
- switch (surrounding_points.size())
- {
- case 2:
- {
- // If two points are passed, these are the two vertices, so
- // we can only compute degree-1 intermediate points.
- for (unsigned int i=0; i<n; ++i)
- {
- const double x = line_support_points.point(i+1)[0];
- w[1] = x;
- w[0] = (1-x);
- points[i] = manifold.get_new_point(surrounding_points, w);
- }
- break;
- }
-
- case 4:
- {
- Assert(spacedim >= 2, ExcImpossibleInDim(spacedim));
- const unsigned m=
- static_cast<unsigned int>(std::sqrt(static_cast<double>(n)));
- // is n a square number
- Assert(m*m==n, ExcInternalError());
-
- // If m points are passed, two of them are the outer vertices, and we
- // can only compute (degree-1)*(degree-1) intermediate points.
- Point<2> p;
- for (unsigned int i=0, c=0; i<m; ++i)
- {
- p[1] = line_support_points.point(1+i)[0];
- for (unsigned int j=0; j<m; ++j, ++c)
- {
- p[0] = line_support_points.point(1+j)[0];
-
- for (unsigned int l=0; l<4; ++l)
- w[l] = GeometryInfo<2>::d_linear_shape_function(p, l);
-
- points[c]=manifold.get_new_point(surrounding_points, w);
- }
- }
- break;
- }
-
- case 8:
- {
- Assert(spacedim >= 3, ExcImpossibleInDim(spacedim));
- unsigned int m=1;
- for ( ; m < n; ++m)
- if (m*m*m == n)
- break;
- // is n a cube number
- Assert(m*m*m==n, ExcInternalError());
-
- Point<3> p;
- for (unsigned int k=0, c=0; k<m; ++k)
- {
- p[2] = line_support_points.point(1+k)[0];
- for (unsigned int i=0; i<m; ++i)
- {
- p[1] = line_support_points.point(1+i)[0];
- for (unsigned int j=0; j<m; ++j, ++c)
- {
- p[0] = line_support_points.point(1+j)[0];
-
- for (unsigned int l=0; l<8; ++l)
- w[l] = GeometryInfo<3>::d_linear_shape_function(p, l);
-
- points[c]=manifold.get_new_point(surrounding_points, w);
- }
- }
- }
- break;
- }
- default:
- Assert(false, ExcInternalError());
- break;
- }
- }
-
-
-
-
/**
* Ask the manifold descriptor to return intermediate points on the object
* pointed to by the TriaIterator @p iter. This function tries to be
* Boundary<dim,spacedim>.
*/
template <int dim, int spacedim, class TriaIterator>
- void get_intermediate_points_on_object(const Manifold<dim, spacedim> &manifold,
- const QGaussLobatto<1> &line_support_points,
- const TriaIterator &iter,
-
std::vector<Point<spacedim> > &points)
+ void add_intermediate_points(const Manifold<dim, spacedim> &manifold,
+ const Table<2,double> &weight_table,
+ const TriaIterator &iter,
+ std::vector<Point<spacedim> > &points)
{
const unsigned int structdim = TriaIterator::AccessorType::structure_dimension;
// Try backward compatibility option.
- if (const Boundary<dim,spacedim> *boundary
- = dynamic_cast<const Boundary<dim,spacedim> *>(&manifold))
+ const Boundary<dim,spacedim> *boundary
+ = dynamic_cast<const Boundary<dim,spacedim> *>(&manifold);
+ if (structdim < spacedim && boundary != NULL)
// This is actually a boundary. Call old methods.
{
+ std::vector<Point<spacedim> > new_points(weight_table.size(0));
switch (structdim)
{
case 1:
{
const typename Triangulation<dim,spacedim>::line_iterator line = iter;
- boundary->get_intermediate_points_on_line(line, points);
- return;
+ boundary->get_intermediate_points_on_line(line, new_points);
+ break;
}
case 2:
{
const typename Triangulation<dim,spacedim>::quad_iterator quad = iter;
- boundary->get_intermediate_points_on_quad(quad, points);
- return;
+ boundary->get_intermediate_points_on_quad(quad, new_points);
+ break;
}
default:
Assert(false, ExcInternalError());
- return;
+ break;
}
+ for (unsigned int i=0; i<new_points.size(); ++i)
+ points.push_back(new_points[i]);
}
else
{
std::vector<Point<spacedim> > sp(GeometryInfo<structdim>::vertices_per_cell);
for (unsigned int i=0; i<sp.size(); ++i)
sp[i] = iter->vertex(i);
- get_intermediate_points(manifold, line_support_points, sp, points);
+ manifold.add_new_points(sp, weight_table, points);
}
}
- /**
- * Take a <tt>support_point_weights_on_hex(quad)</tt> and apply it to the vector
- * @p a to compute the inner support points as a linear combination of the
- * exterior points.
- *
- * The vector @p a initially contains the locations of the @p n_outer
- * points, the @p n_inner computed inner points are appended.
- *
- * See equation (7) of the `mapping' report.
- */
- template <int spacedim>
- void add_weighted_interior_points(const Table<2,double> &lvs,
- std::vector<Point<spacedim> > &a)
+
+ template <typename Iterator>
+ bool check_identical_manifolds_of_quads(const Iterator &)
{
- const unsigned int n_inner_apply=lvs.n_rows();
- const unsigned int n_outer_apply=lvs.n_cols();
- Assert(a.size()==n_outer_apply,
- ExcDimensionMismatch(a.size(), n_outer_apply));
-
- // compute each inner point as linear combination of the outer points. the
- // weights are given by the lvs entries, the outer points are the first
- // (existing) elements of a
- for (unsigned int unit_point=0; unit_point<n_inner_apply; ++unit_point)
- {
- Assert(lvs.n_cols()==n_outer_apply, ExcInternalError());
- Point<spacedim> p;
- for (unsigned int k=0; k<n_outer_apply; ++k)
- p+=lvs[unit_point][k]*a[k];
+ Assert(false, ExcNotImplemented());
+ return true;
+ }
- a.push_back(p);
- }
+ bool check_identical_manifolds_of_quads(const Triangulation<3,3>::cell_iterator &cell)
+ {
+ for (unsigned int f=0; f<GeometryInfo<3>::faces_per_cell; ++f)
+ if (&cell->face(f)->get_manifold() != &cell->get_manifold())
+ return false;
+ return true;
}
}
+
template <int dim, int spacedim>
void
MappingQGeneric<dim,spacedim>::
cell->get_manifold() :
line->get_manifold() );
- get_intermediate_points_on_object (manifold, line_support_points, line, line_points);
+ line_points.resize(0);
+ add_intermediate_points (manifold, support_point_weights_interior[0],
+ line, line_points);
if (dim==3)
{
void
MappingQGeneric<3,3>::
add_quad_support_points(const Triangulation<3,3>::cell_iterator &cell,
- std::vector<Point<3> > &a) const
+ std::vector<Point<3> > &a) const
{
const unsigned int faces_per_cell = GeometryInfo<3>::faces_per_cell;
// used if face quad at boundary or entirely in the interior of the domain
- std::vector<Point<3> > quad_points ((polynomial_degree-1)*(polynomial_degree-1));
+ std::vector<Point<3> > quad_points;
std::vector<Point<3> > tmp_points;
// loop over all faces and collect points on them
// On a quad, we have to check whether the manifold should determine the
// point distribution or rather a weighted sum should be created. This
// is the same logic as in the compute_mapping_support_points function
- // below: if we have a StraightBoundary but no "real" boundary
- // description doing something else, we should go for the smoothing of
- // the support points through add_weighted_interior_points
- if (smooth_support_points_on_flat_manifold == false ||
- std::string(typeid(face->get_manifold()).name()).find("StraightBoundary") ==
- std::string::npos)
+ // below
+ bool all_manifold_ids_are_equal = true;
+ for (unsigned int l=0; l<GeometryInfo<3>::lines_per_face; ++l)
+ if (&face->line(l)->get_manifold() != &face->get_manifold())
+ all_manifold_ids_are_equal = false;
+ if (smooth_support_points == false || all_manifold_ids_are_equal ||
+ (dynamic_cast<const Boundary<3,3> *>(&face->get_manifold()) &&
+ std::string(typeid(face->get_manifold()).name()).find("StraightBoundary") ==
+ std::string::npos))
{
// ask the boundary/manifold object to return intermediate points on it
- get_intermediate_points_on_object(face->get_manifold(), line_support_points,
- face, quad_points);
+ quad_points.resize(0);
+ add_intermediate_points(face->get_manifold(), support_point_weights_interior[1],
+ 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
a[GeometryInfo<3>::vertices_per_cell +
(polynomial_degree-1)*
GeometryInfo<3>::face_to_cell_lines(face_no,line) + i];
- add_weighted_interior_points (support_point_weights_on_quad, tmp_points);
- a.insert(a.end(), tmp_points.begin()+4+4*(polynomial_degree-1),
- tmp_points.end());
+ face->get_manifold().add_new_points (tmp_points,
+ support_point_weights_perimeter_to_interior[1],
+ a);
}
}
}
void
MappingQGeneric<2,3>::
add_quad_support_points(const Triangulation<2,3>::cell_iterator &cell,
- std::vector<Point<3> > &a) const
+ std::vector<Point<3> > &a) const
{
- std::vector<Point<3> > quad_points ((polynomial_degree-1)*(polynomial_degree-1));
- get_intermediate_points_on_object (cell->get_manifold(), line_support_points,
- cell, quad_points);
- for (unsigned int i=0; i<quad_points.size(); ++i)
- a.push_back(quad_points[i]);
+ add_intermediate_points (cell->get_manifold(), support_point_weights_interior[1],
+ cell, a);
}
compute_mapping_support_points(const typename Triangulation<dim,spacedim>::cell_iterator &cell) const
{
// get the vertices first
- std::vector<Point<spacedim> > a(GeometryInfo<dim>::vertices_per_cell);
+ std::vector<Point<spacedim> > a;
+ a.reserve(Utilities::fixed_power<dim>(polynomial_degree+1));
for (unsigned int i=0; i<GeometryInfo<dim>::vertices_per_cell; ++i)
- a[i] = cell->vertex(i);
+ a.push_back(cell->vertex(i));
- if (this->polynomial_degree>1)
- switch (dim)
- {
- case 1:
- add_line_support_points(cell, a);
- break;
- case 2:
- // in 2d, add the points on the four bounding lines to the exterior
- // (outer) points
- add_line_support_points(cell, a);
-
- // then get the support points on the quad. if we have a
- // StraightBoundary but no "real" boundary description doing something
- // else, we should go for the smoothing of the support points through
- // add_weighted_interior_points. we need to go through typeid because
- // curved boundaries are derived from StraightBoundary and thus a
- // plain dynamic_cast will not suffice.
- if (dim != spacedim)
- add_quad_support_points(cell, a);
- else if (smooth_support_points_on_flat_manifold == false ||
- std::string(typeid(cell->get_manifold()).name()).find("StraightBoundary") ==
- std::string::npos)
- {
- std::vector<Point<spacedim> > quad_points (Utilities::fixed_power<dim>(polynomial_degree-1));
- get_intermediate_points_on_object(cell->get_manifold(), line_support_points, cell, quad_points);
- for (unsigned int i=0; i<quad_points.size(); ++i)
- a.push_back(quad_points[i]);
- }
- else
- add_weighted_interior_points (support_point_weights_on_quad, a);
- break;
+ if (this->polynomial_degree > 1)
+ {
+ // check if all entities have the same manifold id which is when we can
+ // simply ask the manifold for all points
+ Assert(dim<=3, ExcImpossibleInDim(dim));
+ bool all_manifold_ids_are_equal = (dim == spacedim);
+ if (dim > 1)
+ for (unsigned int l=0; l<GeometryInfo<dim>::lines_per_cell; ++l)
+ if (&cell->line(l)->get_manifold() != &cell->get_manifold())
+ all_manifold_ids_are_equal = false;
+ if (dim == 3)
+ if (check_identical_manifolds_of_quads(cell) == false)
+ all_manifold_ids_are_equal = false;
+ if (all_manifold_ids_are_equal)
+ {
+ std::vector<Point<spacedim> > vertices(a);
+ cell->get_manifold().add_new_points(vertices, support_point_weights_cell, a);
+ return a;
+ }
- case 3:
- {
- // in 3d also add the points located on the boundary faces
- add_line_support_points (cell, a);
- add_quad_support_points (cell, a);
-
- // then compute the interior points. if we have a StraightBoundary but
- // no "real" boundary description doing something else, we should go
- // for the smoothing of the support points through
- // add_weighted_interior_points. we need to go through typeid because
- // curved boundaries are derived from StraightBoundary and thus a
- // plain dynamic_cast will not suffice.
- if (smooth_support_points_on_flat_manifold == false ||
- std::string(typeid(cell->get_manifold()).name()).find("StraightBoundary") ==
- std::string::npos)
- {
- std::vector<Point<spacedim> > hex_points (Utilities::fixed_power<dim>(polynomial_degree-1));
- get_intermediate_points_on_object(cell->get_manifold(), line_support_points, cell, hex_points);
- for (unsigned int i=0; i<hex_points.size(); ++i)
- a.push_back(hex_points[i]);
- }
- else
- add_weighted_interior_points (support_point_weights_on_hex, a);
- break;
- }
+ switch (dim)
+ {
+ case 1:
+ add_line_support_points(cell, a);
+ break;
+ case 2:
+ // in 2d, add the points on the four bounding lines to the exterior
+ // (outer) points
+ add_line_support_points(cell, a);
+
+ // then get the support points on the quad. decide whether to choose
+ // the smoothed variant that include all points and goes directly
+ // through the manifold or only the vertex points.
+ if (dim != spacedim)
+ add_quad_support_points(cell, a);
+ else if (smooth_support_points == false)
+ add_intermediate_points(cell->get_manifold(), support_point_weights_interior[1],
+ cell, a);
+ else
+ {
+ std::vector<Point<spacedim> > tmp_points(a);
+ cell->get_manifold().add_new_points(tmp_points,
+ support_point_weights_perimeter_to_interior[1],
+ a);
+ }
+ break;
+
+ case 3:
+ // in 3d also add the points located on the boundary faces
+ add_line_support_points (cell, a);
+ add_quad_support_points (cell, a);
+
+ // then compute the interior points. decide whether to choose
+ // the smoothed variant that include all points and goes directly
+ // through the manifold or only the vertex points.
+ if (smooth_support_points == false)
+ add_intermediate_points(cell->get_manifold(), support_point_weights_interior[2],
+ cell, a);
+ else
+ {
+ std::vector<Point<spacedim> > tmp_points(a);
+ cell->get_manifold().add_new_points(tmp_points,
+ support_point_weights_perimeter_to_interior[2],
+ a);
+ }
+ break;
- default:
- Assert(false, ExcNotImplemented());
- break;
- }
+ default:
+ Assert(false, ExcNotImplemented());
+ break;
+ }
+ }
return a;
}
projects / dealii.git / commitdiff