* Virtual destructor.
*/
virtual ~Mapping ();
+
/**
* Prepare internal data
* structures and fill in values
virtual InternalDataBase*
get_face_data (const UpdateFlags flags,
const Quadrature<dim-1>& quadrature) const = 0;
+
/**
* Prepare internal data
* structure for transformation
* @p{p_real} on the real cell
* @p{cell} and returns @p{p_real}.
*/
+//TODO: document meaning of mdata argument
virtual Point<dim> transform_unit_to_real_cell (
const typename Triangulation<dim>::cell_iterator cell,
const Point<dim> &p,
/**
* Constructor.
*/
- InternalData (const Quadrature<dim>& quadrature);
+ InternalData (const Quadrature<dim> &quadrature);
/**
* Length of the cell in
* different coordinate
* directions.
*/
- std::vector<double> length;
+ Tensor<1,dim> length;
/**
* Vector of all quadrature
std::vector<std::vector<Tensor<1,dim> > > unit_tangentials;
/**
- * Auxuliary vectors for internal use.
+ * Auxiliary vectors for internal use.
*/
std::vector<std::vector<Tensor<1,dim> > > aux;
};
InternalData& data,
std::vector<Point<dim> > &quadrature_points,
std::vector<Point<dim> >& normal_vectors) const;
+
+ private:
+ /**
+ * Value to indicate that a given
+ * face or subface number is
+ * invalid.
+ */
+ static const unsigned int invalid_face_number = static_cast<unsigned int>(-1);
};
//TODO: fill_fe_face_values should exist in a version doing all faces
-// to save initialization time.
+//TODO: to save initialization time.
+
+//TODO: doc: laplace_on_quad_vector, compute_laplace_on_quad, etc all
+//TODO: reference to each other, there is no description what they actually
+//TODO: do or are good for
/**
* Mapping class that uses Qp-mappings on boundary AND on inner
*
* For @p{degree==2} this function
* sets the
- * @p{laplace_on_quad_vector} to
+ * @p{laplace_on_hex_vector} to
* the hardcoded data. For
* @p{degree>2} this vector is
* computed.
* function. For @p{dim=1} and 2 this
* function is empty.
*/
+//TODO: rename function to add_quad_support_points, to unify notation
void add_face_support_points(const typename Triangulation<dim>::cell_iterator &cell,
std::vector<Point<dim> > &a) const;
* @p{p_real} on the real cell
* @p{cell} and returns @p{p_real}.
*/
+//TODO: document meaning of mdata argument
virtual Point<dim> transform_unit_to_real_cell (
const typename Triangulation<dim>::cell_iterator cell,
const Point<dim> &p,
get_new_point_on_quad (const typename Triangulation<dim>::quad_iterator &quad) const;
/**
- * Return intermediate points
- * on the boundary line.
- *
+ * Return equally spaced
+ * intermediate points on a
+ * boundary line.
+//TODO: this function is also called for inner lines!?
+ *
* The number of points requested
* is given by the size of the
* vector @p{points}. It is the
typename std::vector<Point<dim> > &points) const;
/**
- * Return intermediate points
- * on the boundary quad.
+ * Return equally spaced
+ * intermediate points on a
+ * boundary quad.
+//TODO: this function is also called for inner quads!?
*
* The number of points requested
* is given by the size of the
* p{n+1} partitions of equal
* lengths.
*
- * Refer to the general documentation of
- * this class and the documentation of the
+ * Refer to the general
+ * documentation of this class
+ * and the documentation of the
* base class.
*/
virtual void
* @p{(m+1)(m+1)} subquads of equal
* size.
*
- * Refer to the general documentation of
- * this class and the documentation of the
+ * Refer to the general
+ * documentation of this class
+ * and the documentation of the
* base class.
*/
virtual void
poly = new TensorProductPolynomials<dim> (v);
}
- // generate permutation/rotation index sets to generate some matrices from others
+ // generate permutation/rotation
+ // index sets to generate some
+ // matrices from others
std::vector<unsigned int> right;
std::vector<unsigned int> top;
rotate_indices (right, 'Z');
// do some internal book-keeping
build_renumbering (*this, degree, renumber);
+//TODO: externalize this to a proper template function
#if deal_II_dimension > 1
build_face_renumbering (FiniteElementData<dim-1>(FE_Q<dim-1>::get_dpo_vector(degree),1),
degree,
#if (deal_II_dimension == 1)
template<>
-const unsigned int Mapping<deal_II_dimension>::normal_directions[2] =
+const unsigned int Mapping<1>::normal_directions[2] =
{
1, 0
};
#if (deal_II_dimension == 2)
template<>
-const unsigned int Mapping<deal_II_dimension>::normal_directions[4] =
+const unsigned int Mapping<2>::normal_directions[4] =
{
2, 0, 3, 1
};
#if (deal_II_dimension == 3)
template<>
-const unsigned int Mapping<deal_II_dimension>::normal_directions[6] =
+const unsigned int Mapping<3>::normal_directions[6] =
{
3, 2, 5, 0, 4, 1
};
{}
-/*------------------------------ InternalData ------------------------------*/
-
-
-template <int dim>
-Mapping<dim>::InternalDataBase::~InternalDataBase ()
-{}
+/*------------------------------ InternalDataBase ------------------------------*/
template <int dim>
+template <int dim>
+Mapping<dim>::InternalDataBase::~InternalDataBase ()
+{}
+
+
/*------------------------------ InternalData ------------------------------*/
+template <int dim>
+const unsigned int MappingCartesian<dim>::invalid_face_number;
+
+
+
template<int dim>
MappingCartesian<dim>::InternalData::InternalData (const Quadrature<dim>& q)
:
- length (dim, 0.),
quadrature_points (q.get_points ())
{}
return update_default;
}
+
+
template <int dim>
UpdateFlags
MappingCartesian<dim>::update_each (const UpdateFlags in) const
}
+
template <int dim>
Mapping<dim>::InternalDataBase*
-MappingCartesian<dim>::get_data (const UpdateFlags flags,
- const Quadrature<dim>& q) const
+MappingCartesian<dim>::get_data (const UpdateFlags flags,
+ const Quadrature<dim> &q) const
{
Assert (flags & update_normal_vectors == 0, ExcNotImplemented());
InternalData* data = new InternalData (q);
template <int dim>
Mapping<dim>::InternalDataBase*
MappingCartesian<dim>::get_face_data (const UpdateFlags,
- const Quadrature<dim-1>& quadrature) const
+ const Quadrature<dim-1>& quadrature) const
{
QProjector<dim> q (quadrature, false);
InternalData* data = new InternalData (q);
template <int dim>
Mapping<dim>::InternalDataBase*
MappingCartesian<dim>::get_subface_data (const UpdateFlags,
- const Quadrature<dim-1>& quadrature) const
+ const Quadrature<dim-1> &quadrature) const
{
QProjector<dim> q (quadrature, true);
InternalData* data = new InternalData (q);
template <int dim>
void
MappingCartesian<dim>::compute_fill (const typename DoFHandler<dim>::cell_iterator &cell,
- const unsigned int face_no,
- const unsigned int sub_no,
- InternalData& data,
+ const unsigned int face_no,
+ const unsigned int sub_no,
+ InternalData &data,
std::vector<Point<dim> > &quadrature_points,
- std::vector<Point<dim> >& normal_vectors) const
+ std::vector<Point<dim> > &normal_vectors) const
{
- const UpdateFlags update_flags(data.current_update_flags());
+//TODO: does it make sense to query the update flags in data if we did not set any in the get_*_data functions?
+ UpdateFlags update_flags(data.current_update_flags());
const unsigned int npts = quadrature_points.size ();
unsigned int offset = 0;
- bool onface = false;
- if (face_no == static_cast<unsigned int> (-1))
+ if (face_no != invalid_face_number)
{
- Assert (sub_no == static_cast<unsigned int> (-1), ExcInternalError());
- } else {
- onface = true;
// Add 1 on both sides of
// assertion to avoid compiler
// warning about testing
// unsigned int < 0 in 1d.
Assert (face_no+1 < GeometryInfo<dim>::faces_per_cell+1,
ExcIndexRange (face_no, 0, GeometryInfo<dim>::faces_per_cell));
- if (sub_no == static_cast<unsigned int> (-1))
+
+ if (sub_no == invalid_face_number)
+ // called from FEFaceValues
offset = face_no * quadrature_points.size();
else
{
+ // called from FESubfaceValue
Assert (sub_no+1 < GeometryInfo<dim>::subfaces_per_face+1,
ExcIndexRange (sub_no, 0, GeometryInfo<dim>::subfaces_per_face));
offset = (face_no * GeometryInfo<dim>::subfaces_per_face + sub_no)
* quadrature_points.size();
}
}
+ else
+ // invalid face number, so
+ // subface should be invalid as
+ // well
+ Assert (sub_no == invalid_face_number, ExcInternalError());
+
+ // let @p{start} be the origin of a
+ // local coordinate system. it is
+ // chosen as the (lower) left
+ // vertex
+ const Point<dim> start = cell->vertex(0);
- // Compute start point and sizes along axes.
- // Strange vertex numbering makes this complicated again.
-
- Point<dim> start = cell->vertex (0);
-
+ // Compute start point and sizes
+ // along axes. Strange vertex
+ // numbering makes this complicated
+ // again.
switch (dim)
{
case 1:
Assert(false, ExcNotImplemented());
}
-
+
+ // transform quadrature point. this
+ // is obtained simply by scaling
+ // unit coordinates with lengths in
+ // each direction
if (update_flags & update_q_points)
{
Assert (quadrature_points.size() == npts,
ExcDimensionMismatch(quadrature_points.size(), npts));
- for (unsigned int i=0;i<npts;++i)
+ for (unsigned int i=0; i<npts; ++i)
{
- Point<dim> p = start;
- for (unsigned int d=0;d<dim;++d)
- p(d) += data.length[d]*data.quadrature_points[i+offset](d);
- quadrature_points[i] = p;
+ quadrature_points[i] = start;
+ for (unsigned int d=0; d<dim; ++d)
+ quadrature_points[i](d) += data.length[d] *
+ data.quadrature_points[i+offset](d);
}
}
-
+
+ // compute normal vectors. since
+ // cells are aligned to coordinate
+ // axes, they are simply vectors
+ // with exactly one entry equal to
+ // 1 or -1
if (update_flags & update_normal_vectors)
{
+ Assert (normal_vectors.size() == npts,
+ ExcDimensionMismatch(normal_vectors.size(), npts));
Point<dim> n;
switch (100*dim+face_no)
{
default:
Assert (false, ExcInternalError());
}
- for (unsigned int i=0;i<npts;++i)
- normal_vectors[i] = n;
+ // furthermore, all normal
+ // vectors on a face are equal
+ std::fill_n (normal_vectors.begin(), normal_vectors.end(), n);
}
data.first_cell = false;
std::vector<Point<dim> >& quadrature_points,
std::vector<double>& JxW_values) const
{
- InternalData *data_ptr = dynamic_cast<InternalData *> (&mapping_data);
- Assert(data_ptr!=0, ExcInternalError());
- InternalData &data=*data_ptr;
+ // convert data object to internal
+ // data for this class. fails with
+ // an exception if that is not
+ // possible
+ InternalData &data = dynamic_cast<InternalData&> (mapping_data);
typename std::vector<Point<dim> > dummy;
- compute_fill (cell, static_cast<unsigned int> (-1), static_cast<unsigned int> (-1),
+ compute_fill (cell, invalid_face_number, invalid_face_number,
data,
quadrature_points,
dummy);
+ // compute Jacobian
+ // determinant. all values are
+ // equal and are the product of the
+ // local lengths in each coordinate
+ // direction
+//TODO: does it make sense to query the update flags in data if we did not set any in the get_*_data functions?
if (data.current_update_flags() & update_JxW_values)
{
double J = data.length[0];
template <int dim>
void
MappingCartesian<dim>::fill_fe_face_values (const typename DoFHandler<dim>::cell_iterator &cell,
- const unsigned int face_no,
- const Quadrature<dim-1> &q,
- typename Mapping<dim>::InternalDataBase &mapping_data,
- std::vector<Point<dim> > &quadrature_points,
- std::vector<double> &JxW_values,
- std::vector<Tensor<1,dim> > &boundary_forms,
- std::vector<Point<dim> > &normal_vectors) const
+ const unsigned int face_no,
+ const Quadrature<dim-1> &q,
+ typename Mapping<dim>::InternalDataBase &mapping_data,
+ std::vector<Point<dim> > &quadrature_points,
+ std::vector<double> &JxW_values,
+ std::vector<Tensor<1,dim> > &boundary_forms,
+ std::vector<Point<dim> > &normal_vectors) const
{
- InternalData *data_ptr = dynamic_cast<InternalData *> (&mapping_data);
- Assert(data_ptr!=0, ExcInternalError());
- InternalData &data=*data_ptr;
+ // convert data object to internal
+ // data for this class. fails with
+ // an exception if that is not
+ // possible
+ InternalData &data = dynamic_cast<InternalData&> (mapping_data);
- compute_fill (cell, face_no, static_cast<unsigned int> (-1),
+ compute_fill (cell, face_no, invalid_face_number,
data,
quadrature_points,
normal_vectors);
+ // first compute Jacobian
+ // determinant, which is simply the
+ // product of the local lengths
+ // since the jacobian is diagonal
double J = 1.;
for (unsigned int d=0;d<dim;++d)
if (d != (normal_directions[face_no]/2))
J *= data.length[d];
+//TODO: does it make sense to query the update flags in data if we did not set any in the get_*_data functions?
if (data.current_update_flags() & update_JxW_values)
{
for (unsigned int i=0; i<JxW_values.size();++i)
JxW_values[i] = J * q.weight(i);
}
+//TODO: does it make sense to query the update flags in data if we did not set any in the get_*_data functions?
if (data.current_update_flags() & update_boundary_forms)
{
for (unsigned int i=0; i<boundary_forms.size();++i)
std::vector<Tensor<1,dim> > &boundary_forms,
std::vector<Point<dim> > &normal_vectors) const
{
- InternalData *data_ptr = dynamic_cast<InternalData *> (&mapping_data);
- Assert(data_ptr!=0, ExcInternalError());
- InternalData &data=*data_ptr;
+ // convert data object to internal
+ // data for this class. fails with
+ // an exception if that is not
+ // possible
+ InternalData &data = dynamic_cast<InternalData&> (mapping_data);
compute_fill (cell, face_no, sub_no,
data,
quadrature_points,
normal_vectors);
+ // first compute Jacobian
+ // determinant, which is simply the
+ // product of the local lengths
+ // since the jacobian is diagonal
double J = 1.;
for (unsigned int d=0;d<dim;++d)
if (d != (normal_directions[face_no]/2))
J *= data.length[d];
+//TODO: does it make sense to query the update flags in data if we did not set any in the get_*_data functions?
if (data.current_update_flags() & update_JxW_values)
{
for (unsigned int i=0; i<JxW_values.size();++i)
JxW_values[i] = J * q.weight(i) / GeometryInfo<dim>::subfaces_per_face;
}
+//TODO: does it make sense to query the update flags in data if we did not set any in the get_*_data functions?
if (data.current_update_flags() & update_boundary_forms)
{
for (unsigned int i=0; i<boundary_forms.size();++i)
template <>
void
MappingCartesian<1>::fill_fe_face_values (const DoFHandler<1>::cell_iterator &,
- const unsigned,
- const Quadrature<0>&,
- Mapping<1>::InternalDataBase&,
- std::vector<Point<1> >&,
- std::vector<double>&,
- std::vector<Tensor<1,1> >&,
- std::vector<Point<1> >&) const
+ const unsigned,
+ const Quadrature<0>&,
+ Mapping<1>::InternalDataBase&,
+ std::vector<Point<1> >&,
+ std::vector<double>&,
+ std::vector<Tensor<1,1> >&,
+ std::vector<Point<1> >&) const
{
Assert(false, ExcNotImplemented());
}
template <>
void
MappingCartesian<1>::fill_fe_subface_values (const DoFHandler<1>::cell_iterator &,
- const unsigned,
- const unsigned,
- const Quadrature<0>&,
- Mapping<1>::InternalDataBase&,
- std::vector<Point<1> >&,
- std::vector<double>&,
- std::vector<Tensor<1,1> >&,
- std::vector<Point<1> >&) const
+ const unsigned,
+ const unsigned,
+ const Quadrature<0>&,
+ Mapping<1>::InternalDataBase&,
+ std::vector<Point<1> >&,
+ std::vector<double>&,
+ std::vector<Tensor<1,1> >&,
+ std::vector<Point<1> >&) const
{
Assert(false, ExcNotImplemented());
}
const Point<dim> &p,
const typename Mapping<dim>::InternalDataBase *const m_data) const
{
+//TODO: wouldn't it be simpler in this function to compute local lengths ourselved, rather than relying on an InternalDataBase object?
+
// If m_data!=0 use this
// InternalData.
//
mdata = dynamic_cast<const InternalData *> (m_data);
Assert(mdata!=0, ExcInternalError());
+//TODO: don't we need to reinit() mdata on the present cell or initialize the mdata->length array in some other way?
+
// use now the InternalData, that
// mdata is pointing to, to compute
- // the point in real space.
+ // the point in real space. it is
+ // simply a scaled version of the
+ // unit cell point shifted by the
+ // (lower) left corner of the cell
Point<dim> p_real = cell->vertex(0);
- for (unsigned int d=0;d<dim;++d)
+ for (unsigned int d=0; d<dim; ++d)
p_real(d) += mdata->length[d]*p(d);
+
+//TODO: memory leak: if initially mdata==0, then an object is created but never deleted. see mapping_q1 for a better implementation
return p_real;
}
const typename Triangulation<dim>::cell_iterator cell,
const Point<dim> &p) const
{
- const Point<dim>& start = cell->vertex (0);
+ const Point<dim> &start = cell->vertex(0);
Point<dim> real = p;
real -= start;
inline
void
MappingCartesian<dim>::contravariant_transformation (std::vector<tensor_> &dst,
- const std::vector<tensor_> &src,
- const Mapping<dim>::InternalDataBase &mapping_data,
- const unsigned int src_offset) const
+ const std::vector<tensor_> &src,
+ const Mapping<dim>::InternalDataBase &mapping_data,
+ const unsigned int src_offset) const
{
Assert(tensor_::dimension==dim && tensor_::rank==1, ExcInvalidData());
- const InternalData* data_ptr = dynamic_cast<const InternalData *> (&mapping_data);
- Assert(data_ptr!=0, ExcInternalError());
- const InternalData &data=*data_ptr;
+
+ // convert data object to internal
+ // data for this class. fails with
+ // an exception if that is not
+ // possible
+ const InternalData &data = dynamic_cast<const InternalData&> (mapping_data);
Assert (data.update_flags & update_contravariant_transformation,
typename FEValuesBase<dim>::ExcAccessToUninitializedField());
- typename std::vector<tensor_>::const_iterator vec = src.begin() + src_offset;
- typename std::vector<tensor_>::iterator result = dst.begin();
- typename std::vector<tensor_>::const_iterator end = dst.end();
+ typename std::vector<tensor_>::const_iterator vec = src.begin() + src_offset;
+ typename std::vector<tensor_>::iterator result = dst.begin();
+ typename std::vector<tensor_>::const_iterator end = dst.end();
+ // simply scale by Jacobian
+ // (which is diagonal here)
while (result!=end)
{
- for (unsigned int d=0;d<dim;++d)
+ for (unsigned int d=0; d<dim; ++d)
(*result)[d] = (*vec)[d]*data.length[d];
- vec++;
- result++;
+ ++vec;
+ ++result;
}
}
+
template <int dim>
template <typename tensor_>
inline
void
MappingCartesian<dim>::covariant_transformation (std::vector<tensor_> &dst,
- const std::vector<tensor_> &src,
- const Mapping<dim>::InternalDataBase &mapping_data,
- const unsigned int src_offset) const
+ const std::vector<tensor_> &src,
+ const Mapping<dim>::InternalDataBase &mapping_data,
+ const unsigned int src_offset) const
{
Assert(tensor_::dimension==dim && tensor_::rank==1, ExcInvalidData());
- const InternalData *data_ptr = dynamic_cast<const InternalData *> (&mapping_data);
- Assert(data_ptr!=0, ExcInternalError());
- const InternalData &data=*data_ptr;
+
+ // convert data object to internal
+ // data for this class. fails with
+ // an exception if that is not
+ // possible
+ const InternalData &data = dynamic_cast<const InternalData&> (mapping_data);
Assert (data.update_flags & update_covariant_transformation,
typename FEValuesBase<dim>::ExcAccessToUninitializedField());
- typename std::vector<tensor_>::const_iterator vec = src.begin() + src_offset;
- typename std::vector<tensor_>::iterator result = dst.begin();
- typename std::vector<tensor_>::const_iterator end = dst.end();
-
+ typename std::vector<tensor_>::const_iterator vec = src.begin() + src_offset;
+ typename std::vector<tensor_>::iterator result = dst.begin();
+ typename std::vector<tensor_>::const_iterator end = dst.end();
+
+ // simply scale by inverse Jacobian
+ // (which is diagonal here)
while (result!=end)
{
for (unsigned int d=0;d<dim;++d)
//----------------------------------------------------------------------//
+// explicit instantiations
template class MappingCartesian<deal_II_dimension>;
template<int dim>
-MappingQ<dim>::InternalData::InternalData (unsigned int n_shape_functions):
+MappingQ<dim>::InternalData::InternalData (const unsigned int n_shape_functions):
MappingQ1<dim>::InternalData(n_shape_functions),
- use_mapping_q1_on_current_cell(false),
- mapping_q1_data(1 << dim)
+//TODO: in 1d, use_mapping_q1_on_current_cell is always false, but should be true.
+ use_mapping_q1_on_current_cell(false),
+ mapping_q1_data(1 << dim)
{
is_mapping_q1_data=false;
}
#if deal_II_dimension == 1
+// in 1d, it is irrelevant which polynomial degree to use, since all
+// cells are scaled linearly
template<>
-MappingQ<1>::MappingQ (unsigned int):
+MappingQ<1>::MappingQ (const unsigned int):
laplace_on_quad_vector(0),
laplace_on_hex_vector(0),
degree(1),
// functions for the Qp mapping of
// cells at the boundary.
std::vector<LagrangeEquidistant> v;
- for (unsigned int i=0;i<=degree;++i)
+ for (unsigned int i=0; i<=degree; ++i)
v.push_back(LagrangeEquidistant(degree,i));
tensor_pols = new TensorProductPolynomials<dim> (v);
renumber.resize(n_shape_functions,0);
std::vector<unsigned int> dpo(dim+1, 1);
for (unsigned int i=1; i<dpo.size(); ++i)
- dpo[i]=dpo[i-1]*(degree-1);
+ dpo[i] = dpo[i-1]*(degree-1);
FiniteElementData<dim> fe_data(dpo, 1);
FE_Q<dim>::build_renumbering (fe_data, p, renumber);
// build laplace_on_quad_vector
if (degree>1)
{
- set_laplace_on_quad_vector(laplace_on_quad_vector);
- if (dim==3)
+ if (dim >= 2)
+ set_laplace_on_quad_vector(laplace_on_quad_vector);
+ if (dim >= 3)
set_laplace_on_hex_vector(laplace_on_hex_vector);
}
}
template<>
void
MappingQ<1>::compute_shapes_virtual (const std::vector<Point<1> > &unit_points,
- MappingQ1<1>::InternalData &data) const
+ MappingQ1<1>::InternalData &data) const
{
MappingQ1<1>::compute_shapes_virtual(unit_points, data);
}
}
+
template <int dim>
Mapping<dim>::InternalDataBase*
MappingQ<dim>::get_data (const UpdateFlags update_flags,
const unsigned int n_original_q_points,
MappingQ1<dim>::InternalData& mapping_q1_data) const
{
- InternalData *data_ptr = dynamic_cast<InternalData *> (&mapping_q1_data);
- Assert(data_ptr!=0, ExcInternalError());
- InternalData &data=*data_ptr;
+ // convert data object to internal
+ // data for this class. fails with
+ // an exception if that is not
+ // possible
+ InternalData &data = dynamic_cast<InternalData&> (mapping_q1_data);
MappingQ1<dim>::compute_face_data(update_flags, q,
n_original_q_points, data);
-
+
+//TODO: externalize this to a proper template function. or: scrap the whole function for 1d since it doesn't make much sense anyway?
#if (deal_II_dimension>1)
if ((data.update_flags & update_normal_vectors)
&& alternative_normals_computation)
void
MappingQ<dim>::fill_fe_values (const DoFHandler<dim>::cell_iterator &cell,
const Quadrature<dim> &q,
- Mapping<dim>::InternalDataBase &mapping_data,
- std::vector<Point<dim> > &quadrature_points,
- std::vector<double> &JxW_values) const
+ Mapping<dim>::InternalDataBase &mapping_data,
+ std::vector<Point<dim> > &quadrature_points,
+ std::vector<double> &JxW_values) const
{
- InternalData *data_ptr = dynamic_cast<InternalData *> (&mapping_data);
- Assert(data_ptr!=0, ExcInternalError());
- InternalData &data=*data_ptr;
+ // convert data object to internal
+ // data for this class. fails with
+ // an exception if that is not
+ // possible
+ InternalData &data = dynamic_cast<InternalData&> (mapping_data);
data.use_mapping_q1_on_current_cell=!(use_mapping_q_on_all_cells
|| cell->has_boundary_lines());
std::vector<Tensor<1,dim> > &exterior_forms,
std::vector<Point<dim> > &normal_vectors) const
{
- InternalData *data_ptr = dynamic_cast<InternalData *> (&mapping_data);
- Assert(data_ptr!=0, ExcInternalError());
- InternalData &data=*data_ptr;
+ // convert data object to internal
+ // data for this class. fails with
+ // an exception if that is not
+ // possible
+ InternalData &data = dynamic_cast<InternalData&> (mapping_data);
+//TODO: shouldn't we ask whether the face is at the boundary, rather than the cell?
data.use_mapping_q1_on_current_cell=!(use_mapping_q_on_all_cells
|| cell->has_boundary_lines());
std::vector<Tensor<1,dim> > &exterior_forms,
std::vector<Point<dim> > &normal_vectors) const
{
- InternalData *data_ptr = dynamic_cast<InternalData *> (&mapping_data);
- Assert(data_ptr!=0, ExcInternalError());
- InternalData &data=*data_ptr;
+ // convert data object to internal
+ // data for this class. fails with
+ // an exception if that is not
+ // possible
+ InternalData &data = dynamic_cast<InternalData&> (mapping_data);
+//TODO: shouldn't we ask whether the face is at the boundary, rather than the cell?
data.use_mapping_q1_on_current_cell=!(use_mapping_q_on_all_cells
|| cell->has_boundary_lines());
Assert(degree>1, ExcInternalError());
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
double const *loqv_ptr=0;
if (degree==2)
{
}
else if (degree==3)
{
- static double loqv3[4*12]
+ static const double loqv3[4*12]
={80/1053., 1/81., 11/1053., 1/81., 25/117., 44/351.,
7/117., 16/351., 7/117., 16/351., 25/117., 44/351.,
1/81., 80/1053., 1/81., 11/1053., 44/351., 25/117.,
loqv_ptr=&loqv3[0];
}
-
+
if (loqv_ptr!=0)
{
+ // precomputed. copy values to
+ // the loqvs array
loqvs.resize(n_inner_2d);
for (unsigned int unit_point=0; unit_point<n_inner_2d; ++unit_point)
{
loqvs[unit_point].resize(n_outer_2d, 0);
- std::vector<double> &loqv=loqvs[unit_point];
for (unsigned int k=0; k<n_outer_2d; ++k)
- loqv[k]=loqv_ptr[unit_point*n_outer_2d+k];
+ loqvs[unit_point][k]=loqv_ptr[unit_point*n_outer_2d+k];
}
}
else
{
+ // not precomputed, then do so now
if (dim==2)
compute_laplace_vector(loqvs);
+//TODO: what if dim==3?
}
-
+ // 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.size(); ++unit_point)
- Assert(fabs(accumulate(
- loqvs[unit_point].begin(),
- loqvs[unit_point].end(),0.)-1)<1e-13,
+ Assert(fabs(std::accumulate(loqvs[unit_point].begin(),
+ loqvs[unit_point].end(),0.)-1)<1e-13,
ExcInternalError());
// TEST output
{
Assert(degree>1, ExcInternalError());
+ // 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
lohvs.resize(n_inner);
for (unsigned int unit_point=0; unit_point<n_inner; ++unit_point)
{
lohvs[unit_point].resize(n_outer, 0);
- std::vector<double> &loqv=lohvs[unit_point];
for (unsigned int k=0; k<n_outer; ++k)
- loqv[k]=lohv_ptr[unit_point*n_outer+k];
+ lohvs[unit_point][k]=lohv_ptr[unit_point*n_outer+k];
}
}
else
+ // 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
for (unsigned int unit_point=0; unit_point<n_inner; ++unit_point)
- Assert(fabs(accumulate(
- lohvs[unit_point].begin(),
- lohvs[unit_point].end(),0.)-1)<1e-13,
+ Assert(fabs(std::accumulate(lohvs[unit_point].begin(),
+ lohvs[unit_point].end(),0.) - 1)<1e-13,
ExcInternalError());
// TEST output
// compute the shape
// gradients at the quadrature
// points on the unit cell
- QGauss4<dim> quadrature;
+ const QGauss4<dim> quadrature;
const unsigned int n_q_points=quadrature.n_quadrature_points;
InternalData quadrature_data(n_shape_functions);
- quadrature_data.shape_derivatives
- .resize(n_shape_functions * n_q_points);
+ quadrature_data.shape_derivatives.resize(n_shape_functions * n_q_points);
compute_shapes(quadrature.get_points(), quadrature_data);
// Compute the stiffness matrix of
MappingQ<dim>::compute_support_points_laplace(const typename Triangulation<dim>::cell_iterator &cell,
std::vector<Point<dim> > &a) const
{
- a.resize(0);
+ a.resize(GeometryInfo<dim>::vertices_per_cell);
// the vertices first
for (unsigned int i=0; i<GeometryInfo<dim>::vertices_per_cell; ++i)
- a.push_back(cell->vertex(i));
+ a[i] = cell->vertex(i);
if (degree>1)
{
std::vector<Point<dim> > &a) const
{
// the vertices first
+// TODO: check for size of 'a' first?
for (unsigned int i=0; i<GeometryInfo<dim>::vertices_per_cell; ++i)
a.push_back(cell->vertex(i));
// (for dim=3)
fill_quad_support_points_simple (cell, a);
- // then the points in cell
- Point<dim> middle;
- compute_midpoint(a, middle);
-
- if (degree==2)
- a.push_back(middle);
- else if (degree==3)
- // The four points in the
- // cell are located at the
- // midpoint between the
- // middle point and the 4
- // vertices
-
- // TODO: better position of
- // points: transform them by
- // a Q2 transformation.
- for (unsigned int i=0; i<GeometryInfo<dim>::vertices_per_cell; ++i)
- a.push_back(middle*2./3.+cell->vertex(vertex_mapping[i])/3.);
- else if (degree==4)
- {
- Assert(a.size()==16, ExcInternalError());
- a.insert(a.end(), 9, Point<dim>());
+ // then the points in cell. for
+ // this we need the midpoint of
+ // the points already in @p{a}
+ const Point<dim> middle = std::accumulate(a.begin(), a.end(), Point<dim>())
+ / a.size();
- const unsigned int inner_map[8]=
- { 0, 1, 2, 5, 8, 7, 6, 3 };
-
-
- // The nine points in the
- // cell are located at the
- // midpoint between the
- // middle point and (the 4
- // vertices and the face
- // midpoints)
-
- a[16+4]=middle;
- for (unsigned int i=0, j=0; i<GeometryInfo<dim>::vertices_per_cell; ++i)
- {
- a[16+inner_map[j++]]=(middle+cell->vertex(i))/2.;
- a[16+inner_map[j++]]=(middle+(cell->vertex(i)+cell->vertex((i+1)%4))/2.)/2.;
- }
- }
- else
- Assert(false, ExcNotImplemented());
- }
+ switch (degree)
+ {
+ case 2:
+ {
+ a.push_back(middle);
+ break;
+ };
+
+ case 3:
+ {
+ // The four points in the
+ // cell are located at
+ // the midpoint between
+ // the middle point and
+ // the 4 vertices
+
+ // TODO: better position of
+ // points: transform them by
+ // a Q2 transformation.
+ for (unsigned int i=0; i<GeometryInfo<dim>::vertices_per_cell; ++i)
+ a.push_back(middle*2./3.+cell->vertex(vertex_mapping[i])/3.);
+ break;
+ };
+
+ case 4:
+ {
+ Assert(a.size()==16, ExcInternalError());
+ a.insert(a.end(), 9, Point<dim>());
+
+ const unsigned int inner_map[8]=
+ { 0, 1, 2, 5, 8, 7, 6, 3 };
+
+
+ // The nine points in the
+ // cell are located at the
+ // midpoint between the
+ // middle point and (the 4
+ // vertices and the face
+ // midpoints)
+
+ a[16+4]=middle;
+ for (unsigned int i=0, j=0; i<GeometryInfo<dim>::vertices_per_cell; ++i)
+ {
+ a[16+inner_map[j++]]=(middle+cell->vertex(i))/2.;
+ a[16+inner_map[j++]]=(middle+(cell->vertex(i)+cell->vertex((i+1)%4))/2.)/2.;
+ }
+ break;
+ };
+
+ default:
+ Assert(false, ExcNotImplemented());
+ };
+ };
+
Assert(a.size()==n_shape_functions, ExcInternalError());
-
-}
-
-
-
-template <int dim> static
-void compute_midpoint(const std::vector<Point<dim> > &points, Point<dim> &m)
-{
- typename std::vector<Point<dim> >::const_iterator p_iter=points.begin();
- const typename std::vector<Point<dim> >::const_iterator p_end=points.end();
- for (; p_iter!=p_end; ++p_iter)
- m+=*p_iter;
- m/=points.size();
}
template <>
void
MappingQ<1>::add_line_support_points (const Triangulation<1>::cell_iterator &,
- std::vector<Point<1> > &) const
-{}
+ std::vector<Point<1> > &) const
+{
+ // there are no points on bounding
+ // lines which are to be added
+}
#endif
template <int dim>
void
MappingQ<dim>::add_line_support_points (const Triangulation<dim>::cell_iterator &cell,
- std::vector<Point<dim> > &a) const
+ std::vector<Point<dim> > &a) const
{
- const Boundary<dim> *boundary;
+ const Boundary<dim> *boundary = 0;
std::vector<Point<dim> > line_points;
if (degree>2)
line_points.resize(degree-1);
-
- Triangulation<dim>::line_iterator line;
+
+ // 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)
{
- line = cell->line(line_no);
+ const typename Triangulation<dim>::line_iterator line = cell->line(line_no);
if (line->at_boundary())
boundary=&line->get_triangulation().get_boundary(line->boundary_indicator());
else
boundary=&straight_boundary;
-
+
+ // if we only need the
+ // midpoint, then ask for
+ // it. otherwise call the more
+ // complicated functions
if (degree==2)
a.push_back(boundary->get_new_point_on_line(line));
else
{
boundary->get_intermediate_points_on_line (line, line_points);
- for (unsigned int i=0; i<line_points.size(); ++i)
- a.push_back(line_points[i]);
+ a.insert (a.end(), line_points.begin(), line_points.end());
}
}
}
#if deal_II_dimension==3
+//TODO: rename function to add_quad_support_points, to unify notation
template<>
void
MappingQ<3>::add_face_support_points(const Triangulation<3>::cell_iterator &cell,
- std::vector<Point<3> > &a) const
+ std::vector<Point<3> > &a) const
{
- const unsigned int faces_per_cell=GeometryInfo<3>::faces_per_cell,
- vertices_per_face=GeometryInfo<3>::vertices_per_face,
- lines_per_face=GeometryInfo<3>::lines_per_face,
- vertices_per_cell=GeometryInfo<3>::vertices_per_cell;
+ const unsigned int faces_per_cell = GeometryInfo<3>::faces_per_cell,
+ vertices_per_face = GeometryInfo<3>::vertices_per_face,
+ lines_per_face = GeometryInfo<3>::lines_per_face,
+ vertices_per_cell = GeometryInfo<3>::vertices_per_cell;
static const unsigned int face_vertex_to_cell_vertex
[faces_per_cell][vertices_per_face]={{0,1,2,3},
{8,7,11,3}};
- Triangulation<3>::face_iterator face;
+ // loop over all faces and collect points on them
for (unsigned int face_no=0; face_no<faces_per_cell; ++face_no)
{
- face=cell->face(face_no);
+ const Triangulation<3>::face_iterator face=cell->face(face_no);
+
for (unsigned int i=0; i<vertices_per_face; ++i)
Assert(face->vertex_index(i)==
cell->vertex_index(face_vertex_to_cell_vertex[face_no][i]),
Assert(face->line(i)==
cell->line(face_line_to_cell_line[face_no][i]),
ExcInternalError());
-
+
+ // if face at boundary, then
+ // ask boundary object to
+ // return intermediate points
+ // on it
if (face->at_boundary())
{
- std::vector<Point<3> > quad_points;
- quad_points.resize((degree-1)*(degree-1));
+ std::vector<Point<3> > quad_points ((degree-1)*(degree-1));
face->get_triangulation().get_boundary(face->boundary_indicator())
.get_intermediate_points_on_quad (face, quad_points);
-
- for (unsigned int i=0; i<quad_points.size(); ++i)
- a.push_back(quad_points[i]);
+ a.insert (a.end(), quad_points.begin(), quad_points.end());
}
else
{
+ // 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())
++lines_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 (lines_at_boundary>0)
{
// sort the points into b
+//TODO: this is not thread-safe!!! b might be used for objects with
+//TODO: different degrees at the same time!
static std::vector<Point<3> > b;
b.resize(4*degree);
Assert(4*degree==vertices_per_face+lines_per_face*(degree-1),
}
else
{
- std::vector<Point<3> > quad_points;
- quad_points.resize((degree-1)*(degree-1));
+ // face is entirely in
+ // the interior. get
+ // intermediate points
+ // from a straight
+ // boundary object
+ std::vector<Point<3> > quad_points ((degree-1)*(degree-1));
straight_boundary.get_intermediate_points_on_quad (face, quad_points);
-
- for (unsigned int i=0; i<quad_points.size(); ++i)
- a.push_back(quad_points[i]);
+ a.insert (a.end(), quad_points.begin(), quad_points.end());
}
}
}
template <>
void
MappingQ<3>::fill_quad_support_points_simple (const Triangulation<3>::cell_iterator &cell,
- std::vector<Point<3> > &a) const
+ std::vector<Point<3> > &a) const
{
- const Boundary<3> *boundary;
+ const Boundary<3> *boundary = 0;
std::vector<Point<3> > quad_points;
Assert(degree>1, ExcInternalError());
quad_points.resize((degree-1)*(degree-1));
- Triangulation<3>::quad_iterator quad;
for (unsigned int quad_no=0; quad_no<GeometryInfo<3>::quads_per_cell; ++quad_no)
{
- quad = cell->face(quad_no);
+ const Triangulation<3>::quad_iterator quad = cell->face(quad_no);
if (quad->at_boundary())
boundary=&quad->get_triangulation().get_boundary(quad->boundary_indicator());
else
boundary=&straight_boundary;
boundary->get_intermediate_points_on_quad (quad, quad_points);
- for (unsigned int i=0; i<quad_points.size(); ++i)
- a.push_back(quad_points[i]);
+ a.insert (a.end(), quad_points.begin(), quad_points.end());
}
}
-
+// explicit instantiation
template class MappingQ<deal_II_dimension>;
mdata=mdata_local;
}
else
+//TODO: can we assume that m_data->mapping_support_points is initialized
mdata = dynamic_cast<const InternalData *> (m_data);
Assert(mdata!=0, ExcInternalError());