DEAL_II_NAMESPACE_OPEN
+
+namespace
+{
+ template <int dim>
+ std::vector<unsigned int>
+ get_dpo_vector (const unsigned int degree)
+ {
+ std::vector<unsigned int> dpo(dim+1, 1U);
+ for (unsigned int i=1; i<dpo.size(); ++i)
+ dpo[i]=dpo[i-1]*(degree-1);
+ return dpo;
+ }
+}
+
namespace internal
{
namespace MappingQ1
return p_unit;
}
- template <int spacedim>
+
+
+ template <int dim, int spacedim>
+ void compute_shape_function_values_general (const unsigned int n_shape_functions,
+ const std::vector<Point<dim> > &unit_points,
+ typename dealii::MappingQGeneric<dim,spacedim>::InternalData &data)
+ {
+ const unsigned int n_points=unit_points.size();
+
+ // Construct the tensor product polynomials used as shape functions for the
+ // Qp mapping of cells at the boundary.
+ const TensorProductPolynomials<dim>
+ tensor_pols (Polynomials::generate_complete_Lagrange_basis(data.line_support_points.get_points()));
+ Assert (n_shape_functions==tensor_pols.n(),
+ ExcInternalError());
+
+ // then also construct the mapping from lexicographic to the Qp shape function numbering
+ const std::vector<unsigned int>
+ renumber (FETools::
+ lexicographic_to_hierarchic_numbering (
+ FiniteElementData<dim> (get_dpo_vector<dim>(data.polynomial_degree), 1,
+ data.polynomial_degree)));
+
+ std::vector<double> values;
+ std::vector<Tensor<1,dim> > grads;
+ if (data.shape_values.size()!=0)
+ {
+ Assert(data.shape_values.size()==n_shape_functions*n_points,
+ ExcInternalError());
+ values.resize(n_shape_functions);
+ }
+ if (data.shape_derivatives.size()!=0)
+ {
+ Assert(data.shape_derivatives.size()==n_shape_functions*n_points,
+ ExcInternalError());
+ grads.resize(n_shape_functions);
+ }
+
+ std::vector<Tensor<2,dim> > grad2;
+ if (data.shape_second_derivatives.size()!=0)
+ {
+ Assert(data.shape_second_derivatives.size()==n_shape_functions*n_points,
+ ExcInternalError());
+ grad2.resize(n_shape_functions);
+ }
+
+ std::vector<Tensor<3,dim> > grad3;
+ if (data.shape_third_derivatives.size()!=0)
+ {
+ Assert(data.shape_third_derivatives.size()==n_shape_functions*n_points,
+ ExcInternalError());
+ grad3.resize(n_shape_functions);
+ }
+
+ std::vector<Tensor<4,dim> > grad4;
+ if (data.shape_fourth_derivatives.size()!=0)
+ {
+ Assert(data.shape_fourth_derivatives.size()==n_shape_functions*n_points,
+ ExcInternalError());
+ grad4.resize(n_shape_functions);
+ }
+
+
+ if (data.shape_values.size()!=0 ||
+ data.shape_derivatives.size()!=0 ||
+ data.shape_second_derivatives.size()!=0 ||
+ data.shape_third_derivatives.size()!=0 ||
+ data.shape_fourth_derivatives.size()!=0 )
+ for (unsigned int point=0; point<n_points; ++point)
+ {
+ tensor_pols.compute(unit_points[point], values, grads, grad2, grad3, grad4);
+
+ if (data.shape_values.size()!=0)
+ for (unsigned int i=0; i<n_shape_functions; ++i)
+ data.shape(point,renumber[i]) = values[i];
+
+ if (data.shape_derivatives.size()!=0)
+ for (unsigned int i=0; i<n_shape_functions; ++i)
+ data.derivative(point,renumber[i]) = grads[i];
+
+ if (data.shape_second_derivatives.size()!=0)
+ for (unsigned int i=0; i<n_shape_functions; ++i)
+ data.second_derivative(point,renumber[i]) = grad2[i];
+
+ if (data.shape_third_derivatives.size()!=0)
+ for (unsigned int i=0; i<n_shape_functions; ++i)
+ data.third_derivative(point,renumber[i]) = grad3[i];
+
+ if (data.shape_fourth_derivatives.size()!=0)
+ for (unsigned int i=0; i<n_shape_functions; ++i)
+ data.fourth_derivative(point,renumber[i]) = grad4[i];
+ }
+ }
+
+
void
- compute_shape_function_values (const unsigned int n_shape_functions,
- const std::vector<Point<1> > &unit_points,
- typename dealii::MappingQGeneric<1,spacedim>::InternalData &data)
+ compute_shape_function_values_hardcode (const unsigned int n_shape_functions,
+ const std::vector<Point<1> > &unit_points,
+ typename dealii::MappingQGeneric<1,1>::InternalData &data)
{
(void)n_shape_functions;
const unsigned int n_points=unit_points.size();
}
if (data.shape_second_derivatives.size()!=0)
{
- // the following may or may not
- // work if dim != spacedim
- Assert (spacedim == 1, ExcNotImplemented());
-
Assert(data.shape_second_derivatives.size()==n_shape_functions*n_points,
ExcInternalError());
data.second_derivative(k,0)[0][0] = 0;
}
if (data.shape_third_derivatives.size()!=0)
{
- // if lower order derivative don't work, neither should this
- Assert (spacedim == 1, ExcNotImplemented());
-
Assert(data.shape_third_derivatives.size()==n_shape_functions*n_points,
ExcInternalError());
}
if (data.shape_fourth_derivatives.size()!=0)
{
- // if lower order derivative don't work, neither should this
- Assert (spacedim == 1, ExcNotImplemented());
-
Assert(data.shape_fourth_derivatives.size()==n_shape_functions*n_points,
ExcInternalError());
}
- template <int spacedim>
void
- compute_shape_function_values (const unsigned int n_shape_functions,
- const std::vector<Point<2> > &unit_points,
- typename dealii::MappingQGeneric<2,spacedim>::InternalData &data)
+ compute_shape_function_values_hardcode (const unsigned int n_shape_functions,
+ const std::vector<Point<2> > &unit_points,
+ typename dealii::MappingQGeneric<2,2>::InternalData &data)
{
+
(void)n_shape_functions;
const unsigned int n_points=unit_points.size();
for (unsigned int k = 0 ; k < n_points ; ++k)
- template <int spacedim>
void
- compute_shape_function_values (const unsigned int n_shape_functions,
- const std::vector<Point<3> > &unit_points,
- typename dealii::MappingQGeneric<3,spacedim>::InternalData &data)
+ compute_shape_function_values_hardcode (const unsigned int n_shape_functions,
+ const std::vector<Point<3> > &unit_points,
+ typename dealii::MappingQGeneric<3,3>::InternalData &data)
{
(void)n_shape_functions;
const unsigned int n_points=unit_points.size();
}
if (data.shape_second_derivatives.size()!=0)
{
- // the following may or may not
- // work if dim != spacedim
- Assert (spacedim == 3, ExcNotImplemented());
-
Assert(data.shape_second_derivatives.size()==n_shape_functions*n_points,
ExcInternalError());
data.second_derivative(k,0)[0][0] = 0;
}
if (data.shape_third_derivatives.size()!=0)
{
- // if lower order derivative don't work, neither should this
- Assert (spacedim == 3, ExcNotImplemented());
-
Assert(data.shape_third_derivatives.size()==n_shape_functions*n_points,
ExcInternalError());
}
if (data.shape_fourth_derivatives.size()!=0)
{
- // if lower order derivative don't work, neither should this
- Assert (spacedim == 3, ExcNotImplemented());
-
Assert(data.shape_fourth_derivatives.size()==n_shape_functions*n_points,
ExcInternalError());
Tensor<4,3> zero;
-namespace
-{
- template <int dim>
- std::vector<unsigned int>
- get_dpo_vector (const unsigned int degree)
- {
- std::vector<unsigned int> dpo(dim+1, 1U);
- for (unsigned int i=1; i<dpo.size(); ++i)
- dpo[i]=dpo[i-1]*(degree-1);
- return dpo;
- }
-}
-
-template<int dim, int spacedim>
+template<>
void
-MappingQGeneric<dim,spacedim>::InternalData::
-compute_shape_function_values (const std::vector<Point<dim> > &unit_points)
+MappingQGeneric<1,1>::InternalData::
+compute_shape_function_values (const std::vector<Point<1> > &unit_points)
{
// if the polynomial degree is one, then we can simplify code a bit
// by using hard-coded shape functions.
- if ((polynomial_degree == 1)
- &&
- (dim == spacedim))
- internal::MappingQ1::compute_shape_function_values<spacedim> (n_shape_functions,
+ if (polynomial_degree == 1)
+ internal::MappingQ1::compute_shape_function_values_hardcode (n_shape_functions,
unit_points, *this);
else
- // otherwise ask an object that describes the polynomial space
{
- const unsigned int n_points=unit_points.size();
-
- // Construct the tensor product polynomials used as shape functions for the
- // Qp mapping of cells at the boundary.
- const TensorProductPolynomials<dim>
- tensor_pols (Polynomials::generate_complete_Lagrange_basis(line_support_points.get_points()));
- Assert (n_shape_functions==tensor_pols.n(),
- ExcInternalError());
-
- // then also construct the mapping from lexicographic to the Qp shape function numbering
- const std::vector<unsigned int>
- renumber (FETools::
- lexicographic_to_hierarchic_numbering (
- FiniteElementData<dim> (get_dpo_vector<dim>(polynomial_degree), 1,
- polynomial_degree)));
-
- std::vector<double> values;
- std::vector<Tensor<1,dim> > grads;
- if (shape_values.size()!=0)
- {
- Assert(shape_values.size()==n_shape_functions*n_points,
- ExcInternalError());
- values.resize(n_shape_functions);
- }
- if (shape_derivatives.size()!=0)
- {
- Assert(shape_derivatives.size()==n_shape_functions*n_points,
- ExcInternalError());
- grads.resize(n_shape_functions);
- }
-
- std::vector<Tensor<2,dim> > grad2;
- if (shape_second_derivatives.size()!=0)
- {
- Assert(shape_second_derivatives.size()==n_shape_functions*n_points,
- ExcInternalError());
- grad2.resize(n_shape_functions);
- }
-
- std::vector<Tensor<3,dim> > grad3;
- if (shape_third_derivatives.size()!=0)
- {
- Assert(shape_third_derivatives.size()==n_shape_functions*n_points,
- ExcInternalError());
- grad3.resize(n_shape_functions);
- }
-
- std::vector<Tensor<4,dim> > grad4;
- if (shape_fourth_derivatives.size()!=0)
- {
- Assert(shape_fourth_derivatives.size()==n_shape_functions*n_points,
- ExcInternalError());
- grad4.resize(n_shape_functions);
- }
-
-
- if (shape_values.size()!=0 ||
- shape_derivatives.size()!=0 ||
- shape_second_derivatives.size()!=0 ||
- shape_third_derivatives.size()!=0 ||
- shape_fourth_derivatives.size()!=0 )
- for (unsigned int point=0; point<n_points; ++point)
- {
- tensor_pols.compute(unit_points[point], values, grads, grad2, grad3, grad4);
-
- if (shape_values.size()!=0)
- for (unsigned int i=0; i<n_shape_functions; ++i)
- shape(point,renumber[i]) = values[i];
-
- if (shape_derivatives.size()!=0)
- for (unsigned int i=0; i<n_shape_functions; ++i)
- derivative(point,renumber[i]) = grads[i];
-
- if (shape_second_derivatives.size()!=0)
- for (unsigned int i=0; i<n_shape_functions; ++i)
- second_derivative(point,renumber[i]) = grad2[i];
+ // otherwise ask an object that describes the polynomial space
+ internal::MappingQ1::compute_shape_function_values_general<1,1>(n_shape_functions,
+ unit_points,*this);
+ }
+}
- if (shape_third_derivatives.size()!=0)
- for (unsigned int i=0; i<n_shape_functions; ++i)
- third_derivative(point,renumber[i]) = grad3[i];
+template<>
+void
+MappingQGeneric<2,2>::InternalData::
+compute_shape_function_values (const std::vector<Point<2> > &unit_points)
+{
+ // if the polynomial degree is one, then we can simplify code a bit
+ // by using hard-coded shape functions.
+ if (polynomial_degree == 1)
+ internal::MappingQ1::compute_shape_function_values_hardcode (n_shape_functions,
+ unit_points, *this);
+ else
+ {
+ // otherwise ask an object that describes the polynomial space
+ internal::MappingQ1::compute_shape_function_values_general<2,2>(n_shape_functions,
+ unit_points,*this);
+ }
+}
- if (shape_fourth_derivatives.size()!=0)
- for (unsigned int i=0; i<n_shape_functions; ++i)
- fourth_derivative(point,renumber[i]) = grad4[i];
- }
+template<>
+void
+MappingQGeneric<3,3>::InternalData::
+compute_shape_function_values (const std::vector<Point<3> > &unit_points)
+{
+ // if the polynomial degree is one, then we can simplify code a bit
+ // by using hard-coded shape functions.
+ if (polynomial_degree == 1)
+ internal::MappingQ1::compute_shape_function_values_hardcode (n_shape_functions,
+ unit_points, *this);
+ else
+ {
+ // otherwise ask an object that describes the polynomial space
+ internal::MappingQ1::compute_shape_function_values_general<3,3>(n_shape_functions,
+ unit_points,*this);
}
}
+template<int dim, int spacedim>
+void
+MappingQGeneric<dim,spacedim>::InternalData::
+compute_shape_function_values (const std::vector<Point<dim> > &unit_points)
+{
+ // for non-matching combinations of dim and spacedim, just run the general
+ // case
+ internal::MappingQ1::compute_shape_function_values_general<dim,spacedim>(n_shape_functions,
+ unit_points,*this);
+}
+
namespace
{