void
compute_manifold_quadrature_weights(const Quadrature<dim> &quadrature);
+ /**
+ * Store vertices internally.
+ */
+ void
+ store_vertices(const typename Triangulation<dim,spacedim>::cell_iterator &cell) const;
+
/**
* Return an estimate (in bytes) or the memory consumption of this object.
*/
virtual std::size_t memory_consumption () const;
/**
- * Store the current cell.
+ * The current cell vertices.
+ *
+ * Computed each.
+ */
+ mutable std::vector<Point<spacedim> > vertices;
+
+ /**
+ * The current cell.
*
* Computed each.
*/
- typename Triangulation<dim,spacedim>::cell_iterator current_cell;
+ mutable typename Triangulation<dim,spacedim>::cell_iterator cell;
/**
* The actual quadrature on the reference cell.
// */
// mutable typename Triangulation<dim,spacedim>::cell_iterator cell_of_current_support_points;
- // /**
- // * The determinant of the Jacobian in each quadrature point. Filled if
- // * #update_volume_elements.
- // */
- // mutable std::vector<double> volume_elements;
+ /**
+ * The determinant of the Jacobian in each quadrature point. Filled if
+ * #update_volume_elements.
+ */
+ mutable std::vector<double> volume_elements;
+
+ /**
+ * A Q1 Finite element, to compute weights.
+ */
+ const FE_Q<dim> fe_q;
};
*/
const FE_Q<dim,spacedim> 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.
- *
- * Sizes: support_point_weights_on_quad.size()= number of inner
- * unit_support_points support_point_weights_on_quad[i].size()= number of
- * outer unit_support_points, i.e. unit_support_points on the boundary of
- * the quad
- *
- * For the definition of this vector see equation (8) of the `mapping'
- * report.
- */
- Table<2,double> support_point_weights_on_quad;
-
- /**
- * 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.
- *
- * For the definition of this vector see equation (8) of the `mapping'
- * report.
- */
- Table<2,double> support_point_weights_on_hex;
-
- /**
- * Return the locations of support points for the mapping. For example, for
- * $Q_1$ mappings these are the vertices, and for higher order polynomial
- * mappings they are the vertices plus interior points on edges, faces, and
- * the cell interior that are placed in consultation with the Manifold
- * description of the domain and its boundary. However, other classes may
- * override this function differently. In particular, the MappingQ1Eulerian
- * class does exactly this by not computing the support points from the
- * geometry of the current cell but instead evaluating an externally given
- * displacement field in addition to the geometry of the cell.
- *
- * The default implementation of this function is appropriate for most
- * cases. It takes the locations of support points on the boundary of the
- * cell from the underlying manifold. Interior support points (ie. support
- * points in quads for 2d, in hexes for 3d) are then 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.
- *
- * The function works its way from the vertices (which it takes from the
- * given cell) via the support points on the line (for which it calls the
- * add_line_support_points() function) and the support points on the quad
- * faces (in 3d, for which it calls the add_quad_support_points() function).
- * It then adds interior support points that are either computed by
- * interpolation from the surrounding points using weights computed by
- * solving a Laplace equation, or if dim<spacedim, it asks the underlying
- * manifold for the locations of interior points.
- */
- // virtual
- // std::vector<Point<spacedim> >
- // compute_mapping_support_points (const typename Triangulation<dim,spacedim>::cell_iterator &cell) const;
-
/**
* Transforms the point @p p on the real cell to the corresponding point on
* the unit cell @p cell by a Newton iteration.
// transform_real_to_unit_cell_internal (const typename Triangulation<dim,spacedim>::cell_iterator &cell,
// const Point<spacedim> &p,
// const Point<dim> &initial_p_unit) const;
-
- /**
- * Make MappingQ a friend since it needs to call the fill_fe_values()
- * functions on its MappingManifold(1) sub-object.
- */
- template <int, int> friend class MappingQ;
-
};
#ifndef DOXYGEN
-// template<int dim, int spacedim>
-// inline
-// const Tensor<1,dim> &
-// MappingManifold<dim,spacedim>::InternalData::derivative (const unsigned int qpoint,
-// const unsigned int shape_nr) const
-// {
-// Assert(qpoint*n_shape_functions + shape_nr < shape_derivatives.size(),
-// ExcIndexRange(qpoint*n_shape_functions + shape_nr, 0,
-// shape_derivatives.size()));
-// return shape_derivatives [qpoint*n_shape_functions + shape_nr];
-// }
+template<int dim, int spacedim>
+inline
+void
+MappingManifold<dim,spacedim>::InternalData::store_vertices (const typename Triangulation<dim,spacedim>::cell_iterator &cell) const
+{
+ vertices.resize(GeometryInfo<dim>::vertices_per_cell);
+ for (unsigned int i=0; i<GeometryInfo<dim>::vertices_per_cell; ++i)
+ vertices[i] = cell->vertex(i);
+ this->cell = cell;
+}
DEAL_II_NAMESPACE_OPEN
template<int dim, int spacedim>
-MappingManifold<dim,spacedim>::InternalData::InternalData ()
+MappingManifold<dim,spacedim>::InternalData::InternalData () :
+ fe_q(1)
{}
// see if we need the (transformation) shape function values
// and/or gradients and resize the necessary arrays
- if (this->update_each & update_quadrature_points)
+ if (this->update_each & (update_quadrature_points | update_contravariant_transformation) )
compute_manifold_quadrature_weights(q);
if (this->update_each & update_covariant_transformation)
template<int dim, int spacedim>
-MappingManifold<dim,spacedim>::MappingManifold ()
- :
+MappingManifold<dim,spacedim>::MappingManifold () :
fe_q(1)
// 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)),
template<int dim, int spacedim>
-MappingManifold<dim,spacedim>::MappingManifold (const MappingManifold<dim,spacedim> &mapping)
- :
+MappingManifold<dim,spacedim>::MappingManifold (const MappingManifold<dim,spacedim> &mapping) :
fe_q(1)
{}
*/
template <int dim, int spacedim>
void
- maybe_compute_q_points (const typename dealii::Triangulation<dim,spacedim>::cell_iterator &cell,
- const typename QProjector<dim>::DataSetDescriptor data_set,
+ maybe_compute_q_points (const typename QProjector<dim>::DataSetDescriptor data_set,
const typename dealii::MappingManifold<dim,spacedim>::InternalData &data,
std::vector<Point<spacedim> > &quadrature_points)
{
const UpdateFlags update_flags = data.update_each;
+ AssertDimension(data.vertices.size(), GeometryInfo<dim>::vertices_per_cell);
if (update_flags & update_quadrature_points)
{
-
AssertDimension(quadrature_points.size(),
data.cell_manifold_quadrature_weights.size());
- std::vector<Point<spacedim> > vertices;
- for (unsigned int v=0; v<GeometryInfo<dim>::vertices_per_cell; ++v)
- vertices.push_back(cell->vertex(v));
-
for (unsigned int point=0; point<quadrature_points.size(); ++point)
{
- quadrature_points[point] = cell->get_manifold().
- get_new_point(Quadrature<spacedim>(vertices,
+ quadrature_points[point] = data.cell->get_manifold().
+ get_new_point(Quadrature<spacedim>(data.vertices,
data.cell_manifold_quadrature_weights[point]));
}
}
*/
template <int dim, int spacedim>
void
- maybe_update_Jacobians (const typename dealii::Triangulation<dim,spacedim>::cell_iterator &cell,
- const typename dealii::QProjector<dim>::DataSetDescriptor data_set,
+ maybe_update_Jacobians (const typename dealii::QProjector<dim>::DataSetDescriptor data_set,
const typename dealii::MappingManifold<dim,spacedim>::InternalData &data)
{
const UpdateFlags update_flags = data.update_each;
std::fill(data.contravariant.begin(), data.contravariant.end(),
DerivativeForm<1,dim,spacedim>());
+ // Cache of weights used to compute points on the reference cell
+ std::vector<double> weights(GeometryInfo<dim>::vertices_per_cell);
+ AssertDimension(GeometryInfo<dim>::vertices_per_cell,
+ data.vertices.size());
for (unsigned int point=0; point<n_q_points; ++point)
{
// Start by figuring out how to compute the direction in
// the reference space:
const Point<dim> &p = data.quad.point(point+data_set);
+ // And get its image on the manifold:
+ const Point<spacedim> P = data.cell->get_manifold().
+ get_new_point(Quadrature<spacedim>(data.vertices,
+ data.cell_manifold_quadrature_weights[point+data_set]));
+
// Always get the maximum length from the point to the
// boundary of the reference element, to compute the
// tangent vectors from the Manifold object
// which is positive, if the coordinate is < .5,
double L = ai > .5 ? -ai: 1-ai;
- data.contravariant[point][i] = cell->get_manifold().get_tangent_vector(p, np)/L;
+ // Get the weights to compute the np point in real space
+ for (unsigned int j=0; j<GeometryInfo<dim>::vertices_per_cell; ++j)
+ weights[j] = data.fe_q.shape_value(j, np);
+
+ Point<spacedim> NP=data.cell->get_manifold().
+ get_new_point(Quadrature<spacedim>(data.vertices, weights));
+
+ Tensor<1,spacedim> T = data.cell->get_manifold().get_tangent_vector(P, NP);
+
+ for (unsigned int d=0; d<spacedim; ++d)
+ data.contravariant[point][i][d] = T[d]/L;
}
}
const unsigned int n_q_points=quadrature.size();
- internal::maybe_compute_q_points<dim,spacedim> (cell,
- QProjector<dim>::DataSetDescriptor::cell (),
+ data.store_vertices(cell);
+
+ internal::maybe_compute_q_points<dim,spacedim> (QProjector<dim>::DataSetDescriptor::cell (),
data,
output_data.quadrature_points);
-// internal::maybe_update_Jacobians<dim,spacedim> (cell_similarity,
-// QProjector<dim>::DataSetDescriptor::cell (),
-// data);
+ internal::maybe_update_Jacobians<dim,spacedim> (QProjector<dim>::DataSetDescriptor::cell (),
+ data);
// const UpdateFlags update_flags = data.update_each;
// const std::vector<double> &weights=quadrature.get_weights();
projects / dealii.git / commitdiff