From: janssen Date: Thu, 9 Aug 2012 11:47:14 +0000 (+0000) Subject: move all boost::any *_nth_derivative() to Tensor X-Git-Url: https://gitweb.dealii.org/cgi-bin/gitweb.cgi?a=commitdiff_plain;h=7800b14da643ac6e9ed517557e5ac984d207181a;p=dealii-svn.git move all boost::any *_nth_derivative() to Tensor *_nth_derivative() and the boost::any functions become internal git-svn-id: https://svn.dealii.org/branches/branch_higher_derivatives@25803 0785d39b-7218-0410-832d-ea1e28bc413d --- diff --git a/deal.II/include/deal.II/base/tensor_product_polynomials.h b/deal.II/include/deal.II/base/tensor_product_polynomials.h index 73c5f21eed..e0c02f596d 100644 --- a/deal.II/include/deal.II/base/tensor_product_polynomials.h +++ b/deal.II/include/deal.II/base/tensor_product_polynomials.h @@ -129,6 +129,7 @@ class TensorProductPolynomials * compute_value(), * compute_grad() or * compute_hessian() + * compute_nth_derivative() * functions, see below, in a * loop over all tensor product * polynomials. @@ -138,6 +139,12 @@ class TensorProductPolynomials std::vector > &grads, std::vector > &hessians) const; + void compute (const Point &unit_point, + std::vector &values, + std::vector > &grads, + std::vector > &hessians, + std::vector > &nth_derivatives) const; + /** * Computes the value of the * ith tensor product @@ -373,6 +380,12 @@ class AnisotropicPolynomials std::vector > &grads, std::vector > &hessians) const; + void compute (const Point &unit_point, + std::vector &values, + std::vector > &grads, + std::vector > &hessians, + std::vector > &nth_derivatives) const; + /** * Computes the value of the * ith tensor product diff --git a/deal.II/include/deal.II/fe/fe.h b/deal.II/include/deal.II/fe/fe.h index 3813fde941..1f43870c3a 100644 --- a/deal.II/include/deal.II/fe/fe.h +++ b/deal.II/include/deal.II/fe/fe.h @@ -637,11 +637,6 @@ class FiniteElement : public Subscriptor, * @param nth_derivative. The return value * @return can only store Tensor. */ - virtual - boost::any - shape_nth_derivative_internal (const unsigned int i, - const Point &p, - const unsigned int nth_derivative) const; template Tensor @@ -663,12 +658,17 @@ class FiniteElement : public Subscriptor, * @param nth_derivative. The return value * @return can only store Tensor. */ - virtual - boost::any + + template + Tensor shape_nth_derivative_component (const unsigned int i, const Point &p, const unsigned int component, - const unsigned int nth_derivative) const; + const unsigned int nth_derivative) const + { + return boost::any_cast > + (shape_nth_derivative_component_internal (i, p, nth_derivative)); + } /** * Check for non-zero values on a @@ -2664,6 +2664,19 @@ class FiniteElement : public Subscriptor, typename Mapping::InternalDataBase &fe_internal, FEValuesData &data) const = 0; + virtual + boost::any + shape_nth_derivative_internal (const unsigned int i, + const Point &p, + const unsigned int nth_derivative) const; + + virtual + boost::any + shape_nth_derivative_component_internal (const unsigned int i, + const Point &p, + const unsigned int component, + const unsigned int nth_derivative) const; + friend class InternalDataBase; friend class FEValuesBase; friend class FEValues; diff --git a/deal.II/include/deal.II/fe/fe_poly.h b/deal.II/include/deal.II/fe/fe_poly.h index 2f95eaa759..238b014e5c 100644 --- a/deal.II/include/deal.II/fe/fe_poly.h +++ b/deal.II/include/deal.II/fe/fe_poly.h @@ -211,11 +211,16 @@ class FE_Poly : public FiniteElement * for more information about the * semantics of this function. */ - virtual - boost::any + + template + Tensor shape_nth_derivative (const unsigned int i, const Point &p, - const unsigned int nth_derivative) const; + const unsigned int nth_derivative) const + { + return boost::any_cast > + (shape_nth_derivative_internal (i, p, nth_derivative)); + } /** * Return the tensor of the @@ -236,12 +241,17 @@ class FE_Poly : public FiniteElement * were called, provided that the * specified component is zero. */ - virtual - boost::any + + template + Tensor shape_nth_derivative_component (const unsigned int i, const Point &p, const unsigned int component, - const unsigned int nth_derivative) const; + const unsigned int nth_derivative) const + { + return boost::any_cast > + (shape_nth_derivative_component_internal (i, p, nth_derivative)); + } protected: virtual @@ -278,6 +288,18 @@ class FE_Poly : public FiniteElement typename Mapping::InternalDataBase &fe_internal, FEValuesData& data) const ; + virtual + boost::any + shape_nth_derivative_internal (const unsigned int i, + const Point &p, + const unsigned int nth_derivative) const; + + virtual + boost::any + shape_nth_derivative_component_internal (const unsigned int i, + const Point &p, + const unsigned int nth_derivative, + const unsigned int component) const; /** * Determine the values that need @@ -424,6 +446,8 @@ class FE_Poly : public FiniteElement std::vector > > shape_7th_derivatives; std::vector > > shape_8th_derivatives; std::vector > > shape_9th_derivatives; + + std::vector > > shape_nth_derivatives; }; /** diff --git a/deal.II/include/deal.II/fe/fe_poly.templates.h b/deal.II/include/deal.II/fe/fe_poly.templates.h index bac8eaf231..cefcbcaa53 100644 --- a/deal.II/include/deal.II/fe/fe_poly.templates.h +++ b/deal.II/include/deal.II/fe/fe_poly.templates.h @@ -114,7 +114,7 @@ FE_Poly::shape_hessian_component (const unsigned int i, template boost::any -FE_Poly::shape_nth_derivative (const unsigned int i, +FE_Poly::shape_nth_derivative_internal (const unsigned int i, const Point &p, const unsigned int nth_derivative) const { @@ -126,7 +126,7 @@ FE_Poly::shape_nth_derivative (const unsigned int i, template boost::any -FE_Poly::shape_nth_derivative_component (const unsigned int i, +FE_Poly::shape_nth_derivative_component_internal (const unsigned int i, const Point &p, const unsigned int component, const unsigned int nth_derivative) const diff --git a/deal.II/include/deal.II/fe/fe_poly_tensor.h b/deal.II/include/deal.II/fe/fe_poly_tensor.h index 14a24361a8..2aed3b6809 100644 --- a/deal.II/include/deal.II/fe/fe_poly_tensor.h +++ b/deal.II/include/deal.II/fe/fe_poly_tensor.h @@ -169,14 +169,22 @@ class FE_PolyTensor : public FiniteElement * vector valued, an exception is * thrown. */ - virtual boost::any shape_nth_derivative (const unsigned int i, + template Tensor shape_nth_derivative (const unsigned int i, const Point &p, - const unsigned int nth_derivative) const; + const unsigned int nth_derivative) const + { + return boost::any_cast > + (shape_nth_derivative_internal (i, p, nth_derivative)); + } - virtual boost::any shape_nth_derivative_component (const unsigned int i, + template Tensor shape_nth_derivative_component (const unsigned int i, const Point &p, const unsigned int component, - const unsigned int nth_derivative) const; + const unsigned int nth_derivative) const + { + return boost::any_cast > + (shape_nth_derivative_component_internal (i, p, nth_derivative)); + } /** * Given flags, @@ -244,6 +252,17 @@ class FE_PolyTensor : public FiniteElement typename Mapping::InternalDataBase &fe_internal, FEValuesData& data) const ; + virtual boost::any shape_nth_derivative_internal (const unsigned int i, + const Point &p, + const unsigned int nth_derivative) const; + + + virtual boost::any shape_nth_derivative_component_internal (const unsigned int i, + const Point &p, + const unsigned int component, + const unsigned int nth_derivative) const; + + /** * Fields of cell-independent * data for FE_PolyTensor. Stores @@ -283,6 +302,16 @@ class FE_PolyTensor : public FiniteElement * point. */ std::vector< std::vector< DerivativeForm<1, dim, spacedim> > > shape_grads; + std::vector > > shape_hessians; + std::vector > > shape_3rd_derivatives; + std::vector > > shape_4th_derivatives; + std::vector > > shape_5th_derivatives; + std::vector > > shape_6th_derivatives; + std::vector > > shape_7th_derivatives; + std::vector > > shape_8th_derivatives; + std::vector > > shape_9th_derivatives; + + std::vector > > shape_nth_derivatives; }; /** diff --git a/deal.II/include/deal.II/fe/fe_values.h b/deal.II/include/deal.II/fe/fe_values.h index 8ab63b280e..69e9bde9d6 100644 --- a/deal.II/include/deal.II/fe/fe_values.h +++ b/deal.II/include/deal.II/fe/fe_values.h @@ -1474,6 +1474,15 @@ class FEValuesData * derivatives. */ typedef std::vector > > HessianVector; + typedef std::vector > > ThirdDerivativeVector; + typedef std::vector > > FourthDerivativeVector; + typedef std::vector > > FifthDerivativeVector; + typedef std::vector > > SixthDerivativeVector; + typedef std::vector > > SeventhDerivativeVector; + typedef std::vector > > EighthDerivativeVector; + typedef std::vector > > NinthDerivativeVector; + + typedef std::vector > > NthDerivativeVector; /** * Store the values of the shape @@ -1503,6 +1512,14 @@ class FEValuesData * 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 @@ -2019,6 +2036,81 @@ class FEValuesBase : protected FEValuesData, const unsigned int point_no, const unsigned int component) const; + /** + * nth_derivative derivatives of + * the function_noth shape function at + * the point_noth 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 @@ -4931,6 +5023,112 @@ FEValuesBase::shape_2nd_derivative_component (const unsigned int i +template +inline +const boost::any & +FEValuesBase::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_derivativeshape_nth_derivatives.size(), + ExcIndexRange (nth_derivative, 0, this->shape_nth_derivatives.size())); + Assert (ishape_nth_derivatives[nth_derivative].size(), + ExcIndexRange (i, 0, this->shape_nth_derivatives[nth_derivative].size())); + Assert (jshape_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 +inline +boost::any +FEValuesBase::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 inline const FiniteElement & diff --git a/deal.II/source/base/tensor_product_polynomials.cc b/deal.II/source/base/tensor_product_polynomials.cc index 7b052f0f67..e0144ccbbb 100644 --- a/deal.II/source/base/tensor_product_polynomials.cc +++ b/deal.II/source/base/tensor_product_polynomials.cc @@ -354,6 +354,149 @@ compute (const Point &p, + +template +void +TensorProductPolynomials:: +compute (const Point &p, + std::vector &values, + std::vector > &grads, + std::vector > &hessians, + std::vector > &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 (that has enough + // fields for any evaluation of values and + // derivatives) + Table<2,Tensor<1,3> > v(dim, polynomials.size()); + { + std::vector tmp (n_values_and_derivatives); + for (unsigned int d=0; d unsigned int TensorProductPolynomials::n() const diff --git a/deal.II/source/fe/fe.cc b/deal.II/source/fe/fe.cc index e8b465211e..f2724c0b6d 100644 --- a/deal.II/source/fe/fe.cc +++ b/deal.II/source/fe/fe.cc @@ -60,7 +60,9 @@ InternalDataBase::initialize_2nd (const FiniteElement *element, // shifted positions std::vector > 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 shift; @@ -306,7 +308,7 @@ FiniteElement::shape_hessian_component (const unsigned int, template boost::any -FiniteElement::shape_nth_derivative (const unsigned int, +FiniteElement::shape_nth_derivative_internal (const unsigned int, const Point &, const unsigned int) const { @@ -318,7 +320,7 @@ FiniteElement::shape_nth_derivative (const unsigned int, template boost::any -FiniteElement::shape_nth_derivative_component(const unsigned int, +FiniteElement::shape_nth_derivative_component_internal(const unsigned int, const Point &, const unsigned int, const unsigned int) const diff --git a/deal.II/source/fe/fe_poly.cc b/deal.II/source/fe/fe_poly.cc index 0601ef0c55..1f50582f01 100644 --- a/deal.II/source/fe/fe_poly.cc +++ b/deal.II/source/fe/fe_poly.cc @@ -51,6 +51,30 @@ FE_Poly,1,2>::fill_fe_values 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) @@ -88,6 +112,30 @@ FE_Poly,2,3>::fill_fe_values 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) @@ -126,6 +174,30 @@ FE_Poly,1,2>::fill_fe_values ( 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) @@ -160,6 +232,30 @@ FE_Poly,2,3>::fill_fe_values 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) diff --git a/deal.II/source/fe/fe_poly_tensor.cc b/deal.II/source/fe/fe_poly_tensor.cc index 63de560b0d..0d428a0b77 100644 --- a/deal.II/source/fe/fe_poly_tensor.cc +++ b/deal.II/source/fe/fe_poly_tensor.cc @@ -238,7 +238,7 @@ FE_PolyTensor::shape_hessian_component (const unsigned int i, template boost::any -FE_PolyTensor::shape_nth_derivative (const unsigned int, const Point &, const unsigned int) const +FE_PolyTensor::shape_nth_derivative_internal (const unsigned int, const Point &, const unsigned int) const { typedef FiniteElement FEE; Assert(false, typename FEE::ExcFENotPrimitive()); @@ -249,7 +249,7 @@ FE_PolyTensor::shape_nth_derivative (const unsigned int, cons template boost::any -FE_PolyTensor::shape_nth_derivative_component (const unsigned int i, +FE_PolyTensor::shape_nth_derivative_component_internal (const unsigned int i, const Point &p, const unsigned int component, const unsigned int nth_derivative) const @@ -307,6 +307,19 @@ FE_PolyTensor::get_data ( std::vector > grads(0); std::vector > hessians(0); + + //lieber jeder vector einzeln?? + std::vector > third_derivatives(0); + std::vector > fourth_derivatives(0); + std::vector > fifth_derivatives(0); + std::vector > sixth_derivatives(0); + std::vector > seventh_derivatives(0); + std::vector > eighth_derivatives(0); + std::vector > ninth_derivatives(0); + + //oder einer für alles?? + std::vector > nth_derivatives(10); + // initialize fields only if really // necessary. otherwise, don't // allocate memory @@ -332,8 +345,86 @@ FE_PolyTensor::get_data ( // 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;idofs_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;idofs_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;idofs_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;idofs_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;idofs_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;idofs_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;idofs_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;idofs_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;idofs_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;idofs_per_cell;++i) + data->shape_9th_derivatives[i].resize (n_q_points); } // Compute shape function values @@ -344,41 +435,61 @@ FE_PolyTensor::get_data ( // functions v_j if (flags & (update_values | update_gradients)) for (unsigned int k=0; kdofs_per_cell; ++i) + data->shape_values[i][k] = values[i]; + else + for (unsigned int i=0; idofs_per_cell; ++i) { - if (inverse_node_matrix.n_cols() == 0) - for (unsigned int i=0; idofs_per_cell; ++i) - data->shape_values[i][k] = values[i]; - else - for (unsigned int i=0; idofs_per_cell; ++i) - { - Tensor<1,dim> add_values; - for (unsigned int j=0; jdofs_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; jdofs_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; idofs_per_cell; ++i) + data->shape_grads[i][k] = grads[i]; + else + for (unsigned int i=0; idofs_per_cell; ++i) { - if (inverse_node_matrix.n_cols() == 0) - for (unsigned int i=0; idofs_per_cell; ++i) - data->shape_grads[i][k] = grads[i]; - else - for (unsigned int i=0; idofs_per_cell; ++i) - { - Tensor<2,dim> add_grads; - for (unsigned int j=0; jdofs_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; jdofs_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; kdofs_per_cell; ++i) + data->shape_values[i][k] = values[i]; + else + for (unsigned int i=0; idofs_per_cell; ++i) + { + Tensor<1,dim> add_values; + for (unsigned int j=0; jdofs_per_cell; ++j) + add_values += inverse_node_matrix(j,i) * values[j]; + data->shape_values[i][k] = add_values; + } + } + return data; } @@ -435,7 +546,6 @@ FE_PolyTensor::fill_fe_values ( for (unsigned int i=0; idofs_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) { diff --git a/deal.II/source/fe/fe_values.cc b/deal.II/source/fe/fe_values.cc index a839d72159..281af89ddb 100644 --- a/deal.II/source/fe/fe_values.cc +++ b/deal.II/source/fe/fe_values.cc @@ -1737,16 +1737,90 @@ FEValuesData::initialize (const unsigned int n_quadrature_p // 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 (n_quadrature_points)); + } if (flags & update_gradients) + { this->shape_gradients.resize (n_nonzero_shape_components, std::vector > (n_quadrature_points)); + this->shape_nth_derivatives.resize(2); + this->shape_nth_derivatives[1].resize (n_nonzero_shape_components, + std::vector (n_quadrature_points)); + } if (flags & update_hessians) + { this->shape_hessians.resize (n_nonzero_shape_components, std::vector > (n_quadrature_points)); + this->shape_nth_derivatives.resize(3); + this->shape_nth_derivatives[2].resize (n_nonzero_shape_components, + std::vector (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 > (n_quadrature_points)); + this->shape_nth_derivatives[3].resize(n_nonzero_shape_components, + std::vector (n_quadrature_points)); + } + + if (flags & update_4th_derivatives) + { + this->shape_4th_derivatives.resize (n_nonzero_shape_components, + std::vector > (n_quadrature_points)); + this->shape_nth_derivatives[4].resize(n_nonzero_shape_components, + std::vector (n_quadrature_points)); + } + + if (flags & update_5th_derivatives) + { + this->shape_5th_derivatives.resize (n_nonzero_shape_components, + std::vector > (n_quadrature_points)); + this->shape_nth_derivatives[5].resize(n_nonzero_shape_components, + std::vector (n_quadrature_points)); + } + + if (flags & update_6th_derivatives) + { + this->shape_6th_derivatives.resize (n_nonzero_shape_components, + std::vector > (n_quadrature_points)); + this->shape_nth_derivatives[6].resize(n_nonzero_shape_components, + std::vector (n_quadrature_points)); + } + + if (flags & update_7th_derivatives) + { + this->shape_7th_derivatives.resize (n_nonzero_shape_components, + std::vector > (n_quadrature_points)); + this->shape_nth_derivatives[7].resize(n_nonzero_shape_components, + std::vector (n_quadrature_points)); + } + + if (flags & update_8th_derivatives) + { + this->shape_8th_derivatives.resize (n_nonzero_shape_components, + std::vector > (n_quadrature_points)); + this->shape_nth_derivatives[8].resize(n_nonzero_shape_components, + std::vector (n_quadrature_points)); + } + + if (flags & update_9th_derivatives) + { + this->shape_9th_derivatives.resize (n_nonzero_shape_components, + std::vector > (n_quadrature_points)); + this->shape_nth_derivatives[9].resize(n_nonzero_shape_components, + std::vector (n_quadrature_points)); + } if (flags & update_quadrature_points) this->quadrature_points.resize(n_quadrature_points);