std::vector<unsigned int>
FE_DGP<dim>::get_dpo_vector(unsigned int deg)
{
- std::vector<unsigned int> dpo(dim+1, static_cast<unsigned int>(0));
+ std::vector<unsigned int> dpo(dim+1, 0U);
dpo[dim] = ++deg;
for (unsigned int i=1;i<dim;++i)
{
std::vector<unsigned int>
FE_DGQ<dim>::get_dpo_vector(unsigned int deg)
{
- std::vector<unsigned int> dpo(dim+1, static_cast<unsigned int>(0));
+ std::vector<unsigned int> dpo(dim+1, 0U);
dpo[dim] = ++deg;
for (unsigned int i=1;i<dim;++i)
dpo[dim] *= deg;
std::vector<unsigned int>
FE_Q<dim>::get_dpo_vector(const unsigned int deg)
{
- std::vector<unsigned int> dpo(dim+1, static_cast<unsigned int>(1));
+ std::vector<unsigned int> dpo(dim+1, 1U);
for (unsigned int i=1; i<dpo.size(); ++i)
dpo[i]=dpo[i-1]*(deg-1);
return dpo;
std::vector<unsigned int>
FE_Q_Hierarchical<dim>::get_dpo_vector(const unsigned int deg)
{
- std::vector<unsigned int> dpo(dim+1, static_cast<unsigned int>(1));
+ std::vector<unsigned int> dpo(dim+1, 1U);
for (unsigned int i=1; i<dpo.size(); ++i)
dpo[i]=dpo[i-1]*(deg-1);
return dpo;
// cells are aligned to coordinate
// axes, they are simply vectors
// with exactly one entry equal to
- // 1 or -1
+ // 1 or -1. Furthermore, all
+ // normals on a face have the same
+ // value
if (update_flags & update_normal_vectors)
{
Assert (normal_vectors.size() == npts,
ExcDimensionMismatch(normal_vectors.size(), npts));
- Point<dim> n;
+ Assert (face_no < GeometryInfo<dim>::faces_per_cell,
+ ExcInternalError());
+
switch (dim)
{
case 2:
{
- switch (face_no)
- {
- case 0:
- n (1) = -1.;
- break;
- case 1:
- n (0) = 1.;
- break;
- case 2:
- n (1) = 1.;
- break;
- case 3:
- n (0) = -1.;
- break;
- default:
- Assert (false, ExcInternalError());
- }
+ static const Point<dim>
+ normals[GeometryInfo<2>::faces_per_cell]
+ = { Point<dim>(0, -1),
+ Point<dim>(1, 0),
+ Point<dim>(0, 1),
+ Point<dim>(-1, 0) };
+ std::fill (normal_vectors.begin(),
+ normal_vectors.end(),
+ normals[face_no]);
break;
}
case 3:
{
- switch (face_no)
- {
- case 0:
- n (1) = -1.;
- break;
- case 1:
- n (1) = 1.;
- break;
- case 2:
- n (2) = -1.;
- break;
- case 3:
- n (0) = 1.;
- break;
- case 4:
- n (2) = 1.;
- break;
- case 5:
- n (0) = -1.;
- break;
- default:
- Assert (false, ExcInternalError());
- }
+ static const Point<dim>
+ normals[GeometryInfo<3>::faces_per_cell]
+ = { Point<dim>(0, -1, 0),
+ Point<dim>(0, 1, 0),
+ Point<dim>(0, 0, -1),
+ Point<dim>(1, 0, 0),
+ Point<dim>(0, 0, 1),
+ Point<dim>(-1, 0, 0) };
+ std::fill (normal_vectors.begin(),
+ normal_vectors.end(),
+ normals[face_no]);
break;
}
default:
Assert (false, ExcNotImplemented());
}
-
- // furthermore, all normal
- // vectors on a face are equal
- std::fill (normal_vectors.begin(), normal_vectors.end(), n);
}
}
template <int dim>
void
-MappingCartesian<dim>::fill_fe_values (const typename DoFHandler<dim>::cell_iterator& cell,
- const Quadrature<dim>& q,
- typename Mapping<dim>::InternalDataBase& mapping_data,
- std::vector<Point<dim> >& quadrature_points,
- std::vector<double>& JxW_values) const
+MappingCartesian<dim>::
+fill_fe_values (const typename DoFHandler<dim>::cell_iterator& cell,
+ const Quadrature<dim>& q,
+ typename Mapping<dim>::InternalDataBase& mapping_data,
+ std::vector<Point<dim> >& quadrature_points,
+ std::vector<double>& JxW_values) const
{
// convert data object to internal
// data for this class. fails with
template<int dim>
void
-MappingQ<dim>::add_quad_support_points(const typename Triangulation<dim>::cell_iterator &,
- std::vector<Point<dim> > &) const
+MappingQ<dim>::
+add_quad_support_points(const typename Triangulation<dim>::cell_iterator &,
+ std::vector<Point<dim> > &) const
{
Assert (dim > 2, ExcImpossibleInDim(dim));
}
template <int dim>
void
-MappingQ<dim>::transform_covariant (Tensor<1,dim> *begin,
- Tensor<1,dim> *end,
- const Tensor<1,dim> *src,
- const typename Mapping<dim>::InternalDataBase &mapping_data) const
+MappingQ<dim>::
+transform_covariant (Tensor<1,dim> *begin,
+ Tensor<1,dim> *end,
+ const Tensor<1,dim> *src,
+ const typename Mapping<dim>::InternalDataBase &mapping_data) const
{
const typename MappingQ1<dim>::InternalData *q1_data =
dynamic_cast<const typename MappingQ1<dim>::InternalData *> (&mapping_data);
else
tensor = data->covariant.begin();
}
+
while (begin!=end)
- {
- contract (*(begin++), *(src++), *(tensor++));
- }
+ contract (*(begin++), *(src++), *(tensor++));
}
else
tensor = data->covariant.begin();
}
+
while (begin!=end)
- {
- contract (*(begin++), *(src++), *(tensor++));
- }
+ contract (*(begin++), *(src++), *(tensor++));
}
template <int dim>
void
-MappingQ<dim>::transform_contravariant (Tensor<1,dim> *begin,
- Tensor<1,dim> *end,
- const Tensor<1,dim> *src,
- const typename Mapping<dim>::InternalDataBase &mapping_data) const
+MappingQ<dim>::
+transform_contravariant (Tensor<1,dim> *begin,
+ Tensor<1,dim> *end,
+ const Tensor<1,dim> *src,
+ const typename Mapping<dim>::InternalDataBase &mapping_data) const
{
const typename MappingQ1<dim>::InternalData *q1_data =
dynamic_cast<const typename MappingQ1<dim>::InternalData *> (&mapping_data);
else
tensor = data->contravariant.begin();
}
+
while (begin!=end)
- {
- contract (*(begin++), *(tensor++), *(src++));
- }
+ contract (*(begin++), *(tensor++), *(src++));
}
else
tensor = data->contravariant.begin();
}
+
while (begin!=end)
- {
- contract (*(begin++), *(tensor++), *(src++));
- }
+ contract (*(begin++), *(tensor++), *(src++));
}
template <int dim>
Point<dim>
-MappingQ<dim>::transform_unit_to_real_cell (
- const typename Triangulation<dim>::cell_iterator &cell,
- const Point<dim> &p) const
+MappingQ<dim>::
+transform_unit_to_real_cell (const typename Triangulation<dim>::cell_iterator &cell,
+ const Point<dim> &p) const
{
// Use the get_data function to
// create an InternalData with data
// transformation shape values
// already computed at point p.
const Quadrature<dim> point_quadrature(p);
- InternalData *mdata=dynamic_cast<InternalData *> (
- get_data(update_transformation_values, point_quadrature));
- Assert(mdata!=0, ExcInternalError());
+ std::auto_ptr<InternalData>
+ mdata (dynamic_cast<InternalData *> (
+ get_data(update_transformation_values, point_quadrature)));
mdata->use_mapping_q1_on_current_cell = !(use_mapping_q_on_all_cells
|| cell->has_boundary_lines());
- typename MappingQ1<dim>::InternalData *p_data=0;
- if (mdata->use_mapping_q1_on_current_cell)
- p_data=&mdata->mapping_q1_data;
- else
- p_data=mdata;
+ typename MappingQ1<dim>::InternalData
+ *p_data = (mdata->use_mapping_q1_on_current_cell ?
+ &mdata->mapping_q1_data :
+ &*mdata);
compute_mapping_support_points(cell, p_data->mapping_support_points);
- const Point<dim> q=this->transform_unit_to_real_cell_internal(*p_data);
- delete mdata;
- return q;
+ return this->transform_unit_to_real_cell_internal(*p_data);
}
template <int dim>
Point<dim>
-MappingQ<dim>::transform_real_to_unit_cell (
- const typename Triangulation<dim>::cell_iterator &cell,
- const Point<dim> &p) const
+MappingQ<dim>::
+transform_real_to_unit_cell (const typename Triangulation<dim>::cell_iterator &cell,
+ const Point<dim> &p) const
{
// first a Newton iteration based
// on a Q1 mapping
- Point<dim> p_unit=MappingQ1<dim>::transform_real_to_unit_cell(cell, p);
+ Point<dim> p_unit = MappingQ1<dim>::transform_real_to_unit_cell(cell, p);
+ // then a Newton iteration based on
+ // the full MappingQ if we need
+ // this
if (cell->has_boundary_lines() || use_mapping_q_on_all_cells)
{
- // then a Newton iteration
- // based on the full MappingQ
const Quadrature<dim> point_quadrature(p_unit);
- InternalData *mdata=dynamic_cast<InternalData *> (
- get_data(update_transformation_values | update_transformation_gradients,
- point_quadrature));
- Assert(mdata!=0, ExcInternalError());
- mdata->use_mapping_q1_on_current_cell=false;
+ std::auto_ptr<InternalData>
+ mdata (dynamic_cast<InternalData *> (
+ get_data(update_transformation_values |
+ update_transformation_gradients,
+ point_quadrature)));
+
+ mdata->use_mapping_q1_on_current_cell = false;
- std::vector<Point<dim> > &points=mdata->mapping_support_points;
- compute_mapping_support_points(cell, points);
+ std::vector<Point<dim> > &points = mdata->mapping_support_points;
+ compute_mapping_support_points (cell, points);
this->transform_real_to_unit_cell_internal(cell, p, *mdata, p_unit);
-
- delete mdata;
}
return p_unit;
template <int dim>
void
-MappingQ1<dim>::transform_contravariant (Tensor<1,dim> *begin,
- Tensor<1,dim> *end,
- const Tensor<1,dim> *src,
- const typename Mapping<dim>::InternalDataBase &mapping_data) const
+MappingQ1<dim>::
+transform_contravariant (Tensor<1,dim> *begin,
+ Tensor<1,dim> *end,
+ const Tensor<1,dim> *src,
+ const typename Mapping<dim>::InternalDataBase &mapping_data) const
{
const InternalData* data_ptr = dynamic_cast<const InternalData *> (&mapping_data);
Assert(data_ptr!=0, ExcInternalError());
tensor = data.contravariant.begin();
while (begin!=end)
- {
- contract (*(begin++), *(tensor++), *(src++));
- }
+ contract (*(begin++), *(tensor++), *(src++));
}
+
template <int dim>
void
-MappingQ1<dim>::transform_contravariant (Tensor<2,dim> *begin,
- Tensor<2,dim> *end,
- const Tensor<2,dim> *src,
- const typename Mapping<dim>::InternalDataBase &mapping_data) const
+MappingQ1<dim>::
+transform_contravariant (Tensor<2,dim> *begin,
+ Tensor<2,dim> *end,
+ const Tensor<2,dim> *src,
+ const typename Mapping<dim>::InternalDataBase &mapping_data) const
{
const InternalData* data_ptr = dynamic_cast<const InternalData *> (&mapping_data);
Assert(data_ptr!=0, ExcInternalError());
tensor = data.contravariant.begin();
while (begin!=end)
- {
- contract (*(begin++), *(tensor++), *(src++));
- }
+ contract (*(begin++), *(tensor++), *(src++));
}
template <int dim>
Point<dim>
-MappingQ1<dim>::transform_unit_to_real_cell (
- const typename Triangulation<dim>::cell_iterator &cell,
- const Point<dim> &p) const
+MappingQ1<dim>::
+transform_unit_to_real_cell (const typename Triangulation<dim>::cell_iterator &cell,
+ const Point<dim> &p) const
{
// Use the get_data function to
// create an InternalData with data
// already computed at point p.
const Quadrature<dim> point_quadrature(p);
- InternalData *mdata=dynamic_cast<InternalData *> (
- get_data(update_transformation_values, point_quadrature));
- Assert(mdata!=0, ExcInternalError());
+ std::auto_ptr<InternalData>
+ mdata (dynamic_cast<InternalData *> (
+ get_data(update_transformation_values, point_quadrature)));
// compute the mapping support
// points
compute_mapping_support_points(cell, mdata->mapping_support_points);
- const Point<dim> p_real=transform_unit_to_real_cell_internal(*mdata);
- delete mdata;
- return p_real;
+ return transform_unit_to_real_cell_internal(*mdata);
}
+
template <int dim>
Point<dim>
-MappingQ1<dim>::transform_unit_to_real_cell_internal (
- const InternalData &data) const
+MappingQ1<dim>::
+transform_unit_to_real_cell_internal (const InternalData &data) const
{
const unsigned int n_mapping_points=data.mapping_support_points.size();
Assert(data.shape_values.size()==n_mapping_points, ExcInternalError());
// compute the point in real space.
Point<dim> p_real;
for (unsigned int i=0; i<data.mapping_support_points.size(); ++i)
- p_real+=data.mapping_support_points[i] * data.shape(0,i);
+ p_real += data.mapping_support_points[i] * data.shape(0,i);
return p_real;
}
+
template <int dim>
Point<dim>
-MappingQ1<dim>::transform_real_to_unit_cell (
- const typename Triangulation<dim>::cell_iterator &cell,
- const Point<dim> &p) const
+MappingQ1<dim>::
+transform_real_to_unit_cell (const typename Triangulation<dim>::cell_iterator &cell,
+ const Point<dim> &p) const
{
// Let the start value of the
// newton iteration be the center
// derivatives already computed at
// point p_unit.
const Quadrature<dim> point_quadrature(p_unit);
- InternalData *mdata=dynamic_cast<InternalData *> (
- MappingQ1<dim>::get_data(update_transformation_values
- | update_transformation_gradients,
- point_quadrature));
- Assert(mdata!=0, ExcInternalError());
-
- MappingQ1<dim>::compute_mapping_support_points(cell, mdata->mapping_support_points);
- Assert(mdata->mapping_support_points.size()==GeometryInfo<dim>::vertices_per_cell,
+ std::auto_ptr<InternalData>
+ mdata(dynamic_cast<InternalData *> (
+ MappingQ1<dim>::get_data(update_transformation_values
+ | update_transformation_gradients,
+ point_quadrature)));
+
+ MappingQ1<dim>::compute_mapping_support_points (cell,
+ mdata->mapping_support_points);
+ Assert(mdata->mapping_support_points.size() ==
+ GeometryInfo<dim>::vertices_per_cell,
ExcInternalError());
// perform the newton iteration.
transform_real_to_unit_cell_internal(cell, p, *mdata, p_unit);
- delete mdata;
return p_unit;
}
template <int dim>
void
-MappingQ1<dim>::transform_real_to_unit_cell_internal (
- const typename Triangulation<dim>::cell_iterator &cell,
- const Point<dim> &p,
- InternalData &mdata,
- Point<dim> &p_unit) const
+MappingQ1<dim>::
+transform_real_to_unit_cell_internal
+(const typename Triangulation<dim>::cell_iterator &cell,
+ const Point<dim> &p,
+ InternalData &mdata,
+ Point<dim> &p_unit) const
{
const unsigned int n_shapes=mdata.shape_values.size();
Assert(n_shapes!=0, ExcInternalError());
compute_shapes(std::vector<Point<dim> > (1, p_unit), mdata);
// f(x)
- p_real=transform_unit_to_real_cell_internal(mdata);
+ p_real = transform_unit_to_real_cell_internal(mdata);
f = p_real-p;
}
}
cell->index(),
shiftmap_dof_handler);
- // We require the cell to be
- // active. This is determined by
- // the user when looping over all
- // active cells in the problem
- // code.
- Assert (dof_cell->active() == true,
- ExcInactiveCell());
+ // We require the cell to be active
+ // since we can only then get nodal
+ // values for the shifts
+ Assert (dof_cell->active() == true, ExcInactiveCell());
// for Q1 elements, the number of
// support points should equal the