* compute_value(),
* compute_grad() or
* compute_hessian()
+ * compute_nth_derivative()
* functions, see below, in a
* loop over all tensor product
* polynomials.
std::vector<Tensor<1,dim> > &grads,
std::vector<Tensor<2,dim> > &hessians) const;
+ void compute (const Point<dim> &unit_point,
+ std::vector<double> &values,
+ std::vector<Tensor<1,dim> > &grads,
+ std::vector<Tensor<2,dim> > &hessians,
+ std::vector<std::vector<boost::any> > &nth_derivatives) const;
+
/**
* Computes the value of the
* <tt>i</tt>th tensor product
std::vector<Tensor<1,dim> > &grads,
std::vector<Tensor<2,dim> > &hessians) const;
+ void compute (const Point<dim> &unit_point,
+ std::vector<double> &values,
+ std::vector<Tensor<1,dim> > &grads,
+ std::vector<Tensor<2,dim> > &hessians,
+ std::vector<std::vector<boost::any> > &nth_derivatives) const;
+
/**
* Computes the value of the
* <tt>i</tt>th tensor product
* @param nth_derivative. The return value
* @return can only store Tensor<nth_derivative,dim>.
*/
- virtual
- boost::any
- shape_nth_derivative_internal (const unsigned int i,
- const Point<dim> &p,
- const unsigned int nth_derivative) const;
template <int n>
Tensor<n,dim>
* @param nth_derivative. The return value
* @return can only store Tensor<nth_derivative,dim>.
*/
- virtual
- boost::any
+
+ template <int n>
+ Tensor<n,dim>
shape_nth_derivative_component (const unsigned int i,
const Point<dim> &p,
const unsigned int component,
- const unsigned int nth_derivative) const;
+ const unsigned int nth_derivative) const
+ {
+ return boost::any_cast<Tensor<n,dim> >
+ (shape_nth_derivative_component_internal (i, p, nth_derivative));
+ }
/**
* Check for non-zero values on a
typename Mapping<dim,spacedim>::InternalDataBase &fe_internal,
FEValuesData<dim,spacedim> &data) const = 0;
+ virtual
+ boost::any
+ shape_nth_derivative_internal (const unsigned int i,
+ const Point<dim> &p,
+ const unsigned int nth_derivative) const;
+
+ virtual
+ boost::any
+ shape_nth_derivative_component_internal (const unsigned int i,
+ const Point<dim> &p,
+ const unsigned int component,
+ const unsigned int nth_derivative) const;
+
friend class InternalDataBase;
friend class FEValuesBase<dim,spacedim>;
friend class FEValues<dim,spacedim>;
* for more information about the
* semantics of this function.
*/
- virtual
- boost::any
+
+ template <int n>
+ Tensor<n,dim>
shape_nth_derivative (const unsigned int i,
const Point<dim> &p,
- const unsigned int nth_derivative) const;
+ const unsigned int nth_derivative) const
+ {
+ return boost::any_cast<Tensor<n,dim> >
+ (shape_nth_derivative_internal (i, p, nth_derivative));
+ }
/**
* Return the tensor of the
* were called, provided that the
* specified component is zero.
*/
- virtual
- boost::any
+
+ template <int n>
+ Tensor<n,dim>
shape_nth_derivative_component (const unsigned int i,
const Point<dim> &p,
const unsigned int component,
- const unsigned int nth_derivative) const;
+ const unsigned int nth_derivative) const
+ {
+ return boost::any_cast<Tensor<n,dim> >
+ (shape_nth_derivative_component_internal (i, p, nth_derivative));
+ }
protected:
virtual
typename Mapping<dim,spacedim>::InternalDataBase &fe_internal,
FEValuesData<dim,spacedim>& data) const ;
+ virtual
+ boost::any
+ shape_nth_derivative_internal (const unsigned int i,
+ const Point<dim> &p,
+ const unsigned int nth_derivative) const;
+
+ virtual
+ boost::any
+ shape_nth_derivative_component_internal (const unsigned int i,
+ const Point<dim> &p,
+ const unsigned int nth_derivative,
+ const unsigned int component) const;
/**
* Determine the values that need
std::vector<std::vector<Tensor<7,dim> > > shape_7th_derivatives;
std::vector<std::vector<Tensor<8,dim> > > shape_8th_derivatives;
std::vector<std::vector<Tensor<9,dim> > > shape_9th_derivatives;
+
+ std::vector<std::vector<std::vector<boost::any> > > shape_nth_derivatives;
};
/**
template <class POLY, int dim, int spacedim>
boost::any
-FE_Poly<POLY,dim,spacedim>::shape_nth_derivative (const unsigned int i,
+FE_Poly<POLY,dim,spacedim>::shape_nth_derivative_internal (const unsigned int i,
const Point<dim> &p,
const unsigned int nth_derivative) const
{
template <class POLY, int dim, int spacedim>
boost::any
-FE_Poly<POLY,dim,spacedim>::shape_nth_derivative_component (const unsigned int i,
+FE_Poly<POLY,dim,spacedim>::shape_nth_derivative_component_internal (const unsigned int i,
const Point<dim> &p,
const unsigned int component,
const unsigned int nth_derivative) const
* vector valued, an exception is
* thrown.
*/
- virtual boost::any shape_nth_derivative (const unsigned int i,
+ template <int n> Tensor<n,dim> shape_nth_derivative (const unsigned int i,
const Point<dim> &p,
- const unsigned int nth_derivative) const;
+ const unsigned int nth_derivative) const
+ {
+ return boost::any_cast<Tensor<n,dim> >
+ (shape_nth_derivative_internal (i, p, nth_derivative));
+ }
- virtual boost::any shape_nth_derivative_component (const unsigned int i,
+ template <int n> Tensor<n,dim> shape_nth_derivative_component (const unsigned int i,
const Point<dim> &p,
const unsigned int component,
- const unsigned int nth_derivative) const;
+ const unsigned int nth_derivative) const
+ {
+ return boost::any_cast<Tensor<n,dim> >
+ (shape_nth_derivative_component_internal (i, p, nth_derivative));
+ }
/**
* Given <tt>flags</tt>,
typename Mapping<dim,spacedim>::InternalDataBase &fe_internal,
FEValuesData<dim,spacedim>& data) const ;
+ virtual boost::any shape_nth_derivative_internal (const unsigned int i,
+ const Point<dim> &p,
+ const unsigned int nth_derivative) const;
+
+
+ virtual boost::any shape_nth_derivative_component_internal (const unsigned int i,
+ const Point<dim> &p,
+ const unsigned int component,
+ const unsigned int nth_derivative) const;
+
+
/**
* Fields of cell-independent
* data for FE_PolyTensor. Stores
* point.
*/
std::vector< std::vector< DerivativeForm<1, dim, spacedim> > > shape_grads;
+ std::vector<std::vector<DerivativeForm<2,dim, spacedim> > > shape_hessians;
+ std::vector<std::vector<DerivativeForm<3,dim, spacedim> > > shape_3rd_derivatives;
+ std::vector<std::vector<DerivativeForm<4,dim, spacedim> > > shape_4th_derivatives;
+ std::vector<std::vector<DerivativeForm<5,dim, spacedim> > > shape_5th_derivatives;
+ std::vector<std::vector<DerivativeForm<6,dim, spacedim> > > shape_6th_derivatives;
+ std::vector<std::vector<DerivativeForm<7,dim, spacedim> > > shape_7th_derivatives;
+ std::vector<std::vector<DerivativeForm<8,dim, spacedim> > > shape_8th_derivatives;
+ std::vector<std::vector<DerivativeForm<9,dim, spacedim> > > shape_9th_derivatives;
+
+ std::vector<std::vector<std::vector<boost::any> > > shape_nth_derivatives;
};
/**
* derivatives.
*/
typedef std::vector<std::vector<Tensor<2,spacedim> > > HessianVector;
+ typedef std::vector<std::vector<Tensor<3,spacedim> > > ThirdDerivativeVector;
+ typedef std::vector<std::vector<Tensor<4,spacedim> > > FourthDerivativeVector;
+ typedef std::vector<std::vector<Tensor<5,spacedim> > > FifthDerivativeVector;
+ typedef std::vector<std::vector<Tensor<6,spacedim> > > SixthDerivativeVector;
+ typedef std::vector<std::vector<Tensor<7,spacedim> > > SeventhDerivativeVector;
+ typedef std::vector<std::vector<Tensor<8,spacedim> > > EighthDerivativeVector;
+ typedef std::vector<std::vector<Tensor<9,spacedim> > > NinthDerivativeVector;
+
+ typedef std::vector<std::vector<std::vector<boost::any> > > NthDerivativeVector;
/**
* Store the values of the shape
* this field.
*/
HessianVector shape_hessians;
+ ThirdDerivativeVector shape_3rd_derivatives;
+ FourthDerivativeVector shape_4th_derivatives;
+ FifthDerivativeVector shape_5th_derivatives;
+ SixthDerivativeVector shape_6th_derivatives;
+ SeventhDerivativeVector shape_7th_derivatives;
+ EighthDerivativeVector shape_8th_derivatives;
+ NinthDerivativeVector shape_9th_derivatives;
+ NthDerivativeVector shape_nth_derivatives;
/**
* Store an array of weights
const unsigned int point_no,
const unsigned int component) const;
+ /**
+ * <tt>nth_derivative</tt> derivatives of
+ * the <tt>function_no</tt>th shape function at
+ * the <tt>point_no</tt>th quadrature point
+ * with respect to real cell
+ * coordinates. If you want to
+ * get the derivatives in one of
+ * the coordinate directions, use
+ * the appropriate function of
+ * the Tensor class to
+ * extract one component. Since
+ * only a reference to the
+ * derivative values is returned,
+ * there should be no major
+ * performance drawback.
+ *
+ * If the shape function is
+ * vector-valued, then this
+ * returns the only non-zero
+ * component. If the shape
+ * function has more than one
+ * non-zero component (i.e. it is
+ * not primitive), then throw an
+ * exception of type
+ * ExcShapeFunctionNotPrimitive. In
+ * that case, use the
+ * shape_nth_derivative_component()
+ * function.
+ *
+ * The same holds for the arguments
+ * of this function as for the
+ * shape_value() function.
+ */
+ const boost::any &
+ shape_nth_derivative (const unsigned int function_no,
+ const unsigned int point_no,
+ const unsigned int nth_derivative) const;
+
+
+ /**
+ * Return one vector component of
+ * the gradient of a shape
+ * function at a quadrature
+ * point. If the finite element
+ * is scalar, then only component
+ * zero is allowed and the return
+ * value equals that of the
+ * shape_nth_derivative()
+ * function. If the finite
+ * element is vector valued but
+ * all shape functions are
+ * primitive (i.e. they are
+ * non-zero in only one
+ * component), then the value
+ * returned by
+ * shape_nth_derivative()
+ * equals that of this function
+ * for exactly one
+ * component. This function is
+ * therefore only of greater
+ * interest if the shape function
+ * is not primitive, but then it
+ * is necessary since the other
+ * function cannot be used.
+ *
+ * The same holds for the arguments
+ * of this function as for the
+ * shape_value_component() function.
+ */
+ boost::any
+ shape_nth_derivative_component (const unsigned int function_no,
+ const unsigned int point_no,
+ const unsigned int component,
+ const unsigned int nth_derivative) const;
+
//@}
/// @name Access to values of global finite element fields
+template <int dim, int spacedim>
+inline
+const boost::any &
+FEValuesBase<dim,spacedim>::shape_nth_derivative (const unsigned int i,
+ const unsigned int j,
+ const unsigned int nth_derivative) const
+{
+ Assert (i < fe->dofs_per_cell,
+ ExcIndexRange (i, 0, fe->dofs_per_cell));
+ Assert (this->update_flags & update_nth_derivatives(nth_derivative),
+ ExcAccessToUninitializedField());
+ Assert (fe->is_primitive (i),
+ ExcShapeFunctionNotPrimitive(i));
+ Assert (nth_derivative<this->shape_nth_derivatives.size(),
+ ExcIndexRange (nth_derivative, 0, this->shape_nth_derivatives.size()));
+ Assert (i<this->shape_nth_derivatives[nth_derivative].size(),
+ ExcIndexRange (i, 0, this->shape_nth_derivatives[nth_derivative].size()));
+ Assert (j<this->shape_nth_derivatives[nth_derivative][0].size(),
+ ExcIndexRange (j, 0, this->shape_nth_derivatives[nth_derivative][0].size()));
+
+ // if the entire FE is primitive,
+ // then we can take a short-cut:
+ if (fe->is_primitive())
+ return this->shape_nth_derivatives[nth_derivative][i][j];
+ else
+ // otherwise, use the mapping
+ // between shape function numbers
+ // and rows. note that by the
+ // assertions above, we know that
+ // this particular shape function
+ // is primitive, so there is no
+ // question to which vector
+ // component the call of this
+ // function refers
+ return this->shape_nth_derivatives[nth_derivative][this->shape_function_to_row_table[i]][j];
+}
+
+
+
+
+template <int dim, int spacedim>
+inline
+boost::any
+FEValuesBase<dim,spacedim>::shape_nth_derivative_component (const unsigned int i,
+ const unsigned int j,
+ const unsigned int component,
+ const unsigned int nth_derivative) const
+{
+ Assert (i < fe->dofs_per_cell,
+ ExcIndexRange (i, 0, fe->dofs_per_cell));
+ Assert (this->update_flags & update_nth_derivatives(nth_derivative),
+ ExcAccessToUninitializedField());
+ Assert (component < fe->n_components(),
+ ExcIndexRange(component, 0, fe->n_components()));
+
+ // if this particular shape
+ // function is primitive, then we
+ // can take a short-cut by checking
+ // whether the requested component
+ // is the only non-zero one (note
+ // that calling
+ // system_to_component_table only
+ // works if the shape function is
+ // primitive):
+ if (fe->is_primitive(i))
+ {
+ if (component == fe->system_to_component_index(i).first)
+ return this->shape_nth_derivatives[nth_derivative][this->shape_function_to_row_table[i]][j];
+ else
+ return boost::any ();
+ }
+ else
+ {
+ // no, this shape function is
+ // not primitive. then we have
+ // to loop over its components
+ // to find the corresponding
+ // row in the arrays of this
+ // object. before that check
+ // whether the shape function
+ // is non-zero at all within
+ // this component:
+ if (fe->get_nonzero_components(i)[component] == false)
+ return boost::any ();
+
+ // count how many non-zero
+ // component the shape function
+ // has before the one we are
+ // looking for, and add this to
+ // the offset of the first
+ // non-zero component of this
+ // shape function in the arrays
+ // we index presently:
+ const unsigned int
+ row = (this->shape_function_to_row_table[i]
+ +
+ std::count (fe->get_nonzero_components(i).begin(),
+ fe->get_nonzero_components(i).begin()+component,
+ true));
+ return this->shape_nth_derivatives[nth_derivative][row][j];
+ }
+}
+
+
+
+
template <int dim, int spacedim>
inline
const FiniteElement<dim,spacedim> &
+
+template <int dim>
+void
+TensorProductPolynomials<dim>::
+compute (const Point<dim> &p,
+ std::vector<double> &values,
+ std::vector<Tensor<1,dim> > &grads,
+ std::vector<Tensor<2,dim> > &hessians,
+ std::vector<std::vector<boost::any> > &nth_derivatives) const
+{
+ Assert (values.size()==n_tensor_pols || values.size()==0,
+ ExcDimensionMismatch2(values.size(), n_tensor_pols, 0));
+ Assert (grads.size()==n_tensor_pols || grads.size()==0,
+ ExcDimensionMismatch2(grads.size(), n_tensor_pols, 0));
+ Assert (hessians.size()==n_tensor_pols|| hessians.size()==0,
+ ExcDimensionMismatch2(hessians.size(), n_tensor_pols, 0));
+
+ for(unsigned int i=0; i<nth_derivatives.size(); ++i)
+ {
+ //TODO: check all the other derivatives BJ
+ //Assert (hessians.size()==n_tensor_pols|| hessians.size()==0,
+ // ExcDimensionMismatch2(hessians.size(), n_tensor_pols, 0));
+ }
+
+ const bool update_values = (values.size() == n_tensor_pols),
+ update_grads = (grads.size()==n_tensor_pols),
+ update_hessians = (hessians.size()==n_tensor_pols),
+ update_3rd_derivatives = (nth_derivatives[3].size()==n_tensor_pols),
+ update_4th_derivatives = (nth_derivatives[4].size()==n_tensor_pols),
+ update_5th_derivatives = (nth_derivatives[5].size()==n_tensor_pols),
+ update_6th_derivatives = (nth_derivatives[6].size()==n_tensor_pols),
+ update_7th_derivatives = (nth_derivatives[7].size()==n_tensor_pols),
+ update_8th_derivatives = (nth_derivatives[8].size()==n_tensor_pols),
+ update_9th_derivatives = (nth_derivatives[9].size()==n_tensor_pols);
+
+ // check how many
+ // values/derivatives we have to
+ // compute
+ unsigned int n_values_and_derivatives = 0;
+ if (update_values)
+ n_values_and_derivatives = 1;
+ if (update_grads)
+ n_values_and_derivatives = 2;
+ if (update_hessians)
+ n_values_and_derivatives = 3;
+ if (update_3rd_derivatives)
+ n_values_and_derivatives = 4;
+ if (update_4th_derivatives)
+ n_values_and_derivatives = 5;
+ if (update_5th_derivatives)
+ n_values_and_derivatives = 6;
+ if (update_6th_derivatives)
+ n_values_and_derivatives = 7;
+ if (update_7th_derivatives)
+ n_values_and_derivatives = 8;
+ if (update_8th_derivatives)
+ n_values_and_derivatives = 9;
+ if (update_9th_derivatives)
+ n_values_and_derivatives = 10;
+
+
+ // compute the values (and derivatives, if
+ // necessary) of all polynomials at this
+ // evaluation point. to avoid many
+ // reallocation, use one std::vector for
+ // polynomial evaluation and store the
+ // result as Tensor<1,3> (that has enough
+ // fields for any evaluation of values and
+ // derivatives)
+ Table<2,Tensor<1,3> > v(dim, polynomials.size());
+ {
+ std::vector<double> tmp (n_values_and_derivatives);
+ for (unsigned int d=0; d<dim; ++d)
+ for (unsigned int i=0; i<polynomials.size(); ++i)
+ {
+ polynomials[i].value(p(d), tmp);
+ for (unsigned int e=0; e<n_values_and_derivatives; ++e)
+ v(d,i)[e] = tmp[e];
+ };
+ }
+
+ for (unsigned int i=0; i<n_tensor_pols; ++i)
+ {
+ // first get the
+ // one-dimensional indices of
+ // this particular tensor
+ // product polynomial
+ unsigned int indices[dim];
+ compute_index (i, indices);
+
+ if (update_values)
+ {
+ values[i] = 1;
+ for (unsigned int x=0; x<dim; ++x)
+ values[i] *= v(x,indices[x])[0];
+ }
+
+ if (update_grads)
+ for (unsigned int d=0; d<dim; ++d)
+ {
+ grads[i][d] = 1.;
+ for (unsigned int x=0; x<dim; ++x)
+ grads[i][d] *= v(x,indices[x])[d==x];
+ }
+
+ if (update_hessians)
+ for (unsigned int d1=0; d1<dim; ++d1)
+ for (unsigned int d2=0; d2<dim; ++d2)
+ {
+ hessians[i][d1][d2] = 1.;
+ for (unsigned int x=0; x<dim; ++x)
+ {
+ unsigned int derivative=0;
+ if (d1==x || d2==x)
+ {
+ if (d1==d2)
+ derivative=2;
+ else
+ derivative=1;
+ }
+ hessians[i][d1][d2]
+ *= v(x,indices[x])[derivative];
+ }
+ }
+ if (update_3rd_derivatives)
+ {/*do something clever here*/}
+ if (update_4th_derivatives)
+ {/*do something clever here*/}
+ if (update_5th_derivatives)
+ {/*do something clever here*/}
+ if (update_6th_derivatives)
+ {/*do something clever here*/}
+ if (update_7th_derivatives)
+ {/*do something clever here*/}
+ if (update_8th_derivatives)
+ {/*do something clever here*/}
+ if (update_9th_derivatives)
+ {/*do something clever here*/}
+ }
+}
+
+
+
template <int dim>
unsigned int
TensorProductPolynomials<dim>::n() const
// shifted positions
std::vector<Point<dim> > diff_points (quadrature.size());
+ deallog << "before resize" << std::endl;
differences.resize(2*dim);
+ deallog << "after resize " << differences.size() << std::endl;
for (unsigned int d=0; d<dim; ++d)
{
Point<dim> shift;
template <int dim, int spacedim>
boost::any
-FiniteElement<dim,spacedim>::shape_nth_derivative (const unsigned int,
+FiniteElement<dim,spacedim>::shape_nth_derivative_internal (const unsigned int,
const Point<dim> &,
const unsigned int) const
{
template <int dim, int spacedim>
boost::any
-FiniteElement<dim,spacedim>::shape_nth_derivative_component(const unsigned int,
+FiniteElement<dim,spacedim>::shape_nth_derivative_component_internal(const unsigned int,
const Point<dim> &,
const unsigned int,
const unsigned int) const
if (flags & update_gradients && cell_similarity != CellSimilarity::translation)
mapping.transform(fe_data.shape_gradients[k], data.shape_gradients[k],
mapping_data, mapping_covariant);
+
+ if (flags & update_hessians && cell_similarity != CellSimilarity::translation)
+ { /*do something clever*/ }
+
+ if (flags & update_3rd_derivatives && cell_similarity != CellSimilarity::translation)
+ { /*do something clever*/ }
+
+ if (flags & update_4th_derivatives && cell_similarity != CellSimilarity::translation)
+ { /*do something clever*/ }
+
+ if (flags & update_5th_derivatives && cell_similarity != CellSimilarity::translation)
+ { /*do something clever*/ }
+
+ if (flags & update_6th_derivatives && cell_similarity != CellSimilarity::translation)
+ { /*do something clever*/ }
+
+ if (flags & update_7th_derivatives && cell_similarity != CellSimilarity::translation)
+ { /*do something clever*/ }
+
+ if (flags & update_8th_derivatives && cell_similarity != CellSimilarity::translation)
+ { /*do something clever*/ }
+
+ if (flags & update_9th_derivatives && cell_similarity != CellSimilarity::translation)
+ { /*do something clever*/ }
}
if (flags & update_hessians && cell_similarity != CellSimilarity::translation)
if (flags & update_gradients && cell_similarity != CellSimilarity::translation)
mapping.transform(fe_data.shape_gradients[k], data.shape_gradients[k],
mapping_data, mapping_covariant);
+
+ if (flags & update_hessians && cell_similarity != CellSimilarity::translation)
+ { /*do something clever*/ }
+
+ if (flags & update_3rd_derivatives && cell_similarity != CellSimilarity::translation)
+ { /*do something clever*/ }
+
+ if (flags & update_4th_derivatives && cell_similarity != CellSimilarity::translation)
+ { /*do something clever*/ }
+
+ if (flags & update_5th_derivatives && cell_similarity != CellSimilarity::translation)
+ { /*do something clever*/ }
+
+ if (flags & update_6th_derivatives && cell_similarity != CellSimilarity::translation)
+ { /*do something clever*/ }
+
+ if (flags & update_7th_derivatives && cell_similarity != CellSimilarity::translation)
+ { /*do something clever*/ }
+
+ if (flags & update_8th_derivatives && cell_similarity != CellSimilarity::translation)
+ { /*do something clever*/ }
+
+ if (flags & update_9th_derivatives && cell_similarity != CellSimilarity::translation)
+ { /*do something clever*/ }
}
if (flags & update_hessians && cell_similarity != CellSimilarity::translation)
if (flags & update_gradients && cell_similarity != CellSimilarity::translation)
mapping.transform(fe_data.shape_gradients[k], data.shape_gradients[k],
mapping_data, mapping_covariant);
+
+ if (flags & update_hessians && cell_similarity != CellSimilarity::translation)
+ { /*do something clever*/ }
+
+ if (flags & update_3rd_derivatives && cell_similarity != CellSimilarity::translation)
+ { /*do something clever*/ }
+
+ if (flags & update_4th_derivatives && cell_similarity != CellSimilarity::translation)
+ { /*do something clever*/ }
+
+ if (flags & update_5th_derivatives && cell_similarity != CellSimilarity::translation)
+ { /*do something clever*/ }
+
+ if (flags & update_6th_derivatives && cell_similarity != CellSimilarity::translation)
+ { /*do something clever*/ }
+
+ if (flags & update_7th_derivatives && cell_similarity != CellSimilarity::translation)
+ { /*do something clever*/ }
+
+ if (flags & update_8th_derivatives && cell_similarity != CellSimilarity::translation)
+ { /*do something clever*/ }
+
+ if (flags & update_9th_derivatives && cell_similarity != CellSimilarity::translation)
+ { /*do something clever*/ }
}
if (flags & update_hessians && cell_similarity != CellSimilarity::translation)
if (flags & update_gradients && cell_similarity != CellSimilarity::translation)
mapping.transform(fe_data.shape_gradients[k], data.shape_gradients[k],
mapping_data, mapping_covariant);
+
+ if (flags & update_hessians && cell_similarity != CellSimilarity::translation)
+ { /*do something clever*/ }
+
+ if (flags & update_3rd_derivatives && cell_similarity != CellSimilarity::translation)
+ { /*do something clever*/ }
+
+ if (flags & update_4th_derivatives && cell_similarity != CellSimilarity::translation)
+ { /*do something clever*/ }
+
+ if (flags & update_5th_derivatives && cell_similarity != CellSimilarity::translation)
+ { /*do something clever*/ }
+
+ if (flags & update_6th_derivatives && cell_similarity != CellSimilarity::translation)
+ { /*do something clever*/ }
+
+ if (flags & update_7th_derivatives && cell_similarity != CellSimilarity::translation)
+ { /*do something clever*/ }
+
+ if (flags & update_8th_derivatives && cell_similarity != CellSimilarity::translation)
+ { /*do something clever*/ }
+
+ if (flags & update_9th_derivatives && cell_similarity != CellSimilarity::translation)
+ { /*do something clever*/ }
}
if (flags & update_hessians && cell_similarity != CellSimilarity::translation)
template <class POLY, int dim, int spacedim>
boost::any
-FE_PolyTensor<POLY,dim,spacedim>::shape_nth_derivative (const unsigned int, const Point<dim> &, const unsigned int) const
+FE_PolyTensor<POLY,dim,spacedim>::shape_nth_derivative_internal (const unsigned int, const Point<dim> &, const unsigned int) const
{
typedef FiniteElement<dim,spacedim> FEE;
Assert(false, typename FEE::ExcFENotPrimitive());
template <class POLY, int dim, int spacedim>
boost::any
-FE_PolyTensor<POLY,dim,spacedim>::shape_nth_derivative_component (const unsigned int i,
+FE_PolyTensor<POLY,dim,spacedim>::shape_nth_derivative_component_internal (const unsigned int i,
const Point<dim> &p,
const unsigned int component,
const unsigned int nth_derivative) const
std::vector<Tensor<2,dim> > grads(0);
std::vector<Tensor<3,dim> > hessians(0);
+
+ //lieber jeder vector einzeln??
+ std::vector<Tensor<4,dim> > third_derivatives(0);
+ std::vector<Tensor<5,dim> > fourth_derivatives(0);
+ std::vector<Tensor<6,dim> > fifth_derivatives(0);
+ std::vector<Tensor<7,dim> > sixth_derivatives(0);
+ std::vector<Tensor<8,dim> > seventh_derivatives(0);
+ std::vector<Tensor<9,dim> > eighth_derivatives(0);
+ std::vector<Tensor<10,dim> > ninth_derivatives(0);
+
+ //oder einer für alles??
+ std::vector<std::vector<boost::any> > nth_derivatives(10);
+
// initialize fields only if really
// necessary. otherwise, don't
// allocate memory
// that
if (flags & update_hessians)
{
-// hessians.resize (this->dofs_per_cell);
+ hessians.resize (this->dofs_per_cell);
+ //old version BJ
data->initialize_2nd (this, mapping, quadrature);
+
+ //new version BJ
+ //data->shape_hessians.resize (this->dofs_per_cell);
+ //for (unsigned int i=0;i<this->dofs_per_cell;++i)
+ // data->shape_hessians[i].resize (n_q_points);
+ }
+
+ if (flags & update_3rd_derivatives)
+ {
+ nth_derivatives[3].resize (this->dofs_per_cell);
+ data->shape_3rd_derivatives.resize (this->dofs_per_cell);
+ for (unsigned int i=0;i<this->dofs_per_cell;++i)
+ data->shape_3rd_derivatives[i].resize (n_q_points);
+ }
+
+ if (flags & update_4th_derivatives)
+ {
+ nth_derivatives[4].resize (this->dofs_per_cell);
+ data->shape_4th_derivatives.resize (this->dofs_per_cell);
+ for (unsigned int i=0;i<this->dofs_per_cell;++i)
+ data->shape_4th_derivatives[i].resize (n_q_points);
+ }
+
+ if (flags & update_5th_derivatives)
+ {
+ nth_derivatives[5].resize (this->dofs_per_cell);
+ data->shape_5th_derivatives.resize (this->dofs_per_cell);
+ for (unsigned int i=0;i<this->dofs_per_cell;++i)
+ data->shape_5th_derivatives[i].resize (n_q_points);
+ }
+
+ if (flags & update_4th_derivatives)
+ {
+ nth_derivatives[4].resize (this->dofs_per_cell);
+ data->shape_4th_derivatives.resize (this->dofs_per_cell);
+ for (unsigned int i=0;i<this->dofs_per_cell;++i)
+ data->shape_4th_derivatives[i].resize (n_q_points);
+ }
+
+ if (flags & update_5th_derivatives)
+ {
+ nth_derivatives[5].resize (this->dofs_per_cell);
+ data->shape_5th_derivatives.resize (this->dofs_per_cell);
+ for (unsigned int i=0;i<this->dofs_per_cell;++i)
+ data->shape_5th_derivatives[i].resize (n_q_points);
+ }
+
+ if (flags & update_6th_derivatives)
+ {
+ nth_derivatives[6].resize (this->dofs_per_cell);
+ data->shape_6th_derivatives.resize (this->dofs_per_cell);
+ for (unsigned int i=0;i<this->dofs_per_cell;++i)
+ data->shape_6th_derivatives[i].resize (n_q_points);
+ }
+
+ if (flags & update_7th_derivatives)
+ {
+ nth_derivatives[7].resize (this->dofs_per_cell);
+ data->shape_7th_derivatives.resize (this->dofs_per_cell);
+ for (unsigned int i=0;i<this->dofs_per_cell;++i)
+ data->shape_7th_derivatives[i].resize (n_q_points);
+ }
+
+ if (flags & update_8th_derivatives)
+ {
+ nth_derivatives[8].resize (this->dofs_per_cell);
+ data->shape_8th_derivatives.resize (this->dofs_per_cell);
+ for (unsigned int i=0;i<this->dofs_per_cell;++i)
+ data->shape_8th_derivatives[i].resize (n_q_points);
+ }
+
+ if (flags & update_9th_derivatives)
+ {
+ nth_derivatives[9].resize (this->dofs_per_cell);
+ data->shape_9th_derivatives.resize (this->dofs_per_cell);
+ for (unsigned int i=0;i<this->dofs_per_cell;++i)
+ data->shape_9th_derivatives[i].resize (n_q_points);
}
// Compute shape function values
// functions v_j
if (flags & (update_values | update_gradients))
for (unsigned int k=0; k<n_q_points; ++k)
- {
- poly_space.compute(quadrature.point(k),
- values, grads, hessians);
+ {
+ poly_space.compute(quadrature.point(k),
+ values, grads, hessians);
- if (flags & update_values)
+ if (flags & update_values)
+ {
+ if (inverse_node_matrix.n_cols() == 0)
+ for (unsigned int i=0; i<this->dofs_per_cell; ++i)
+ data->shape_values[i][k] = values[i];
+ else
+ for (unsigned int i=0; i<this->dofs_per_cell; ++i)
{
- if (inverse_node_matrix.n_cols() == 0)
- for (unsigned int i=0; i<this->dofs_per_cell; ++i)
- data->shape_values[i][k] = values[i];
- else
- for (unsigned int i=0; i<this->dofs_per_cell; ++i)
- {
- Tensor<1,dim> add_values;
- for (unsigned int j=0; j<this->dofs_per_cell; ++j)
- add_values += inverse_node_matrix(j,i) * values[j];
- data->shape_values[i][k] = add_values;
- }
+ Tensor<1,dim> add_values;
+ for (unsigned int j=0; j<this->dofs_per_cell; ++j)
+ add_values += inverse_node_matrix(j,i) * values[j];
+ data->shape_values[i][k] = add_values;
}
+ }
- if (flags & update_gradients)
+ if (flags & update_gradients)
+ {
+ if (inverse_node_matrix.n_cols() == 0)
+ for (unsigned int i=0; i<this->dofs_per_cell; ++i)
+ data->shape_grads[i][k] = grads[i];
+ else
+ for (unsigned int i=0; i<this->dofs_per_cell; ++i)
{
- if (inverse_node_matrix.n_cols() == 0)
- for (unsigned int i=0; i<this->dofs_per_cell; ++i)
- data->shape_grads[i][k] = grads[i];
- else
- for (unsigned int i=0; i<this->dofs_per_cell; ++i)
- {
- Tensor<2,dim> add_grads;
- for (unsigned int j=0; j<this->dofs_per_cell; ++j)
- add_grads += inverse_node_matrix(j,i) * grads[j];
- data->shape_grads[i][k] = add_grads;
- }
+ Tensor<2,dim> add_grads;
+ for (unsigned int j=0; j<this->dofs_per_cell; ++j)
+ add_grads += inverse_node_matrix(j,i) * grads[j];
+ data->shape_grads[i][k] = add_grads;
}
-
}
+ }
+
+ //new version BJ
+ if (flags & update_hessians)
+ for (unsigned int k=0; k<n_q_points; ++k)
+ {
+ poly_space.compute(quadrature.point(k),
+ values, grads, hessians);
+
+ if (inverse_node_matrix.n_cols() == 0)
+ for (unsigned int i=0; i<this->dofs_per_cell; ++i)
+ data->shape_values[i][k] = values[i];
+ else
+ for (unsigned int i=0; i<this->dofs_per_cell; ++i)
+ {
+ Tensor<1,dim> add_values;
+ for (unsigned int j=0; j<this->dofs_per_cell; ++j)
+ add_values += inverse_node_matrix(j,i) * values[j];
+ data->shape_values[i][k] = add_values;
+ }
+ }
+
return data;
}
for (unsigned int i=0; i<this->dofs_per_cell; ++i)
{
const unsigned int first = data.shape_function_to_row_table[i];
-
if (flags & update_values && cell_similarity != CellSimilarity::translation)
switch (mapping_type)
{
// that we will need to their
// correct size
if (flags & update_values)
+ {
this->shape_values.reinit(n_nonzero_shape_components,
n_quadrature_points);
+ this->shape_nth_derivatives.resize(1);
+ this->shape_nth_derivatives[0].resize(n_nonzero_shape_components,
+ std::vector<boost::any> (n_quadrature_points));
+ }
if (flags & update_gradients)
+ {
this->shape_gradients.resize (n_nonzero_shape_components,
std::vector<Tensor<1,spacedim> > (n_quadrature_points));
+ this->shape_nth_derivatives.resize(2);
+ this->shape_nth_derivatives[1].resize (n_nonzero_shape_components,
+ std::vector<boost::any> (n_quadrature_points));
+ }
if (flags & update_hessians)
+ {
this->shape_hessians.resize (n_nonzero_shape_components,
std::vector<Tensor<2,spacedim> > (n_quadrature_points));
+ this->shape_nth_derivatives.resize(3);
+ this->shape_nth_derivatives[2].resize (n_nonzero_shape_components,
+ std::vector<boost::any> (n_quadrature_points));
+ }
+
+ if (flags & update_derivatives(3,9))
+ this->shape_nth_derivatives.resize(10);
+
+ if (flags & update_3rd_derivatives)
+ {
+ this->shape_3rd_derivatives.resize (n_nonzero_shape_components,
+ std::vector<Tensor<3,spacedim> > (n_quadrature_points));
+ this->shape_nth_derivatives[3].resize(n_nonzero_shape_components,
+ std::vector<boost::any> (n_quadrature_points));
+ }
+
+ if (flags & update_4th_derivatives)
+ {
+ this->shape_4th_derivatives.resize (n_nonzero_shape_components,
+ std::vector<Tensor<4,spacedim> > (n_quadrature_points));
+ this->shape_nth_derivatives[4].resize(n_nonzero_shape_components,
+ std::vector<boost::any> (n_quadrature_points));
+ }
+
+ if (flags & update_5th_derivatives)
+ {
+ this->shape_5th_derivatives.resize (n_nonzero_shape_components,
+ std::vector<Tensor<5,spacedim> > (n_quadrature_points));
+ this->shape_nth_derivatives[5].resize(n_nonzero_shape_components,
+ std::vector<boost::any> (n_quadrature_points));
+ }
+
+ if (flags & update_6th_derivatives)
+ {
+ this->shape_6th_derivatives.resize (n_nonzero_shape_components,
+ std::vector<Tensor<6,spacedim> > (n_quadrature_points));
+ this->shape_nth_derivatives[6].resize(n_nonzero_shape_components,
+ std::vector<boost::any> (n_quadrature_points));
+ }
+
+ if (flags & update_7th_derivatives)
+ {
+ this->shape_7th_derivatives.resize (n_nonzero_shape_components,
+ std::vector<Tensor<7,spacedim> > (n_quadrature_points));
+ this->shape_nth_derivatives[7].resize(n_nonzero_shape_components,
+ std::vector<boost::any> (n_quadrature_points));
+ }
+
+ if (flags & update_8th_derivatives)
+ {
+ this->shape_8th_derivatives.resize (n_nonzero_shape_components,
+ std::vector<Tensor<8,spacedim> > (n_quadrature_points));
+ this->shape_nth_derivatives[8].resize(n_nonzero_shape_components,
+ std::vector<boost::any> (n_quadrature_points));
+ }
+
+ if (flags & update_9th_derivatives)
+ {
+ this->shape_9th_derivatives.resize (n_nonzero_shape_components,
+ std::vector<Tensor<9,spacedim> > (n_quadrature_points));
+ this->shape_nth_derivatives[9].resize(n_nonzero_shape_components,
+ std::vector<boost::any> (n_quadrature_points));
+ }
if (flags & update_quadrature_points)
this->quadrature_points.resize(n_quadrature_points);