From: wolf Date: Sun, 8 Nov 1998 20:36:10 +0000 (+0000) Subject: Work on the implementation of second derivatives. X-Git-Url: https://gitweb.dealii.org/cgi-bin/gitweb.cgi?a=commitdiff_plain;h=a5b9d559ac4e38749f407a3711f5856fac192758;p=dealii-svn.git Work on the implementation of second derivatives. git-svn-id: https://svn.dealii.org/trunk@657 0785d39b-7218-0410-832d-ea1e28bc413d --- diff --git a/deal.II/deal.II/include/fe/fe.h b/deal.II/deal.II/include/fe/fe.h index 69a3bc511d..836f4c892a 100644 --- a/deal.II/deal.II/include/fe/fe.h +++ b/deal.II/deal.II/include/fe/fe.h @@ -589,12 +589,15 @@ class FiniteElement : public FiniteElementBase { * and the given quadrature points on the * unit cell. The Jacobian matrix is to * be computed at every quadrature point. + * The derivative of the jacobian matrix + * is the derivative with respect to the + * unit cell coordinates. * This function has to be in the finite * element class, since different finite * elements need different transformations * of the unit cell to a real cell. * - * The computation of the three fields may + * The computation of these fields may * share some common code, which is why we * put it in one function. However, it may * not always be necessary to really @@ -607,25 +610,33 @@ class FiniteElement : public FiniteElementBase { * of the Jacobi matrix and of the various * structures to be filled. * - * It is provided for the finite element - * class in one space dimension, but for - * higher dimensions, it depends on the - * present fe and needs reimplementation - * by the user. This is due to the fact + * This function is provided for + * the finite element class in + * one space dimension, but for + * higher dimensions, it depends + * on the present fe and needs + * reimplementation by the + * user. This is due to the fact * that the user may want to use - * iso- or subparametric mappings of the - * unit cell to the real cell, which - * makes things much more complicated. + * iso- or subparametric mappings + * of the unit cell to the real + * cell, which makes things much + * more complicated. * - * The #shape_values/grads_transform# - * arrays store the values and gradients - * of the transformation basis functions. - * While this information is not necessary - * for the computation of the other fields, - * it allows for significant speedups, since - * the values and gradients of the transform - * functions at the quadrature points - * need not be recomputed each time this + * The + * #shape_values/grads_transform# + * arrays store the values and + * gradients of the + * transformation basis + * functions. While this + * information is not necessary + * for the computation of the + * other fields, it allows for + * significant speedups, since + * the values and gradients of + * the transform functions at the + * quadrature points need not be + * recomputed each time this * function is called. * * The function assumes that the fields @@ -635,17 +646,18 @@ class FiniteElement : public FiniteElementBase { * This function is more or less an * interface to the #FEValues# class and * should not be used by users unless - * absolutely needed. - */ + * absolutely needed. */ virtual void fill_fe_values (const DoFHandler::cell_iterator &cell, - const vector > &unit_points, - vector > &jacobians, - const bool compute_jacobians, - vector > &support_points, - const bool compute_support_points, - vector > &q_points, - const bool compute_q_points, - const dFMatrix &shape_values_transform, + const vector > &unit_points, + vector > &jacobians, + const bool compute_jacobians, + vector > &jacobians_grad, + const bool compute_jacobians_grad, + vector > &support_points, + const bool compute_support_points, + vector > &q_points, + const bool compute_q_points, + const dFMatrix &shape_values_transform, const vector > > &shape_grads_transform, const Boundary &boundary) const; @@ -758,7 +770,9 @@ class FiniteElement : public FiniteElementBase { const vector > &unit_points, const vector > &global_unit_points, vector > &jacobians, - const bool compute_jacobians, + const bool compute_jacobians, + vector > &jacobians_grad, + const bool compute_jacobians_grad, vector > &support_points, const bool compute_support_points, vector > &q_points, @@ -805,7 +819,9 @@ class FiniteElement : public FiniteElementBase { const vector > &unit_points, const vector > &global_unit_points, vector > &jacobians, - const bool compute_jacobians, + const bool compute_jacobians, + vector > &jacobians_grad, + const bool compute_jacobians_grad, vector > &q_points, const bool compute_q_points, vector &face_jacobi_determinants, diff --git a/deal.II/deal.II/include/fe/fe_lib.criss_cross.h b/deal.II/deal.II/include/fe/fe_lib.criss_cross.h index 68d41b8ea3..6e7d2c4931 100644 --- a/deal.II/deal.II/include/fe/fe_lib.criss_cross.h +++ b/deal.II/deal.II/include/fe/fe_lib.criss_cross.h @@ -333,12 +333,14 @@ class FECrissCross : public FiniteElement { virtual void fill_fe_values (const DoFHandler::cell_iterator &cell, const vector > &unit_points, vector > &jacobians, - const bool compute_jacobians, - vector > &support_points, - const bool compute_support_points, - vector > &q_points, - const bool compute_q_points, - const dFMatrix &shape_values_transform, + const bool compute_jacobians, + vector > &jacobians_grad, + const bool compute_jacobians_grad, + vector > &support_points, + const bool compute_support_points, + vector > &q_points, + const bool compute_q_points, + const dFMatrix &shape_values_transform, const vector > > &shape_grad_transform, const Boundary &boundary) const; diff --git a/deal.II/deal.II/include/fe/fe_linear_mapping.h b/deal.II/deal.II/include/fe/fe_linear_mapping.h index 9be39b199f..0b0af49d04 100644 --- a/deal.II/deal.II/include/fe/fe_linear_mapping.h +++ b/deal.II/deal.II/include/fe/fe_linear_mapping.h @@ -136,12 +136,14 @@ class FELinearMapping : public FiniteElement { virtual void fill_fe_values (const DoFHandler::cell_iterator &cell, const vector > &unit_points, vector > &jacobians, - const bool compute_jacobians, - vector > &support_points, - const bool compute_support_points, - vector > &q_points, - const bool compute_q_points, - const dFMatrix &shape_values_transform, + const bool compute_jacobians, + vector > &jacobians_grad, + const bool compute_jacobians_grad, + vector > &support_points, + const bool compute_support_points, + vector > &q_points, + const bool compute_q_points, + const dFMatrix &shape_values_transform, const vector > > &shape_grad_transform, const Boundary &boundary) const; }; diff --git a/deal.II/deal.II/include/fe/fe_update_flags.h b/deal.II/deal.II/include/fe/fe_update_flags.h index 9e61b96271..793211f699 100644 --- a/deal.II/deal.II/include/fe/fe_update_flags.h +++ b/deal.II/deal.II/include/fe/fe_update_flags.h @@ -71,7 +71,7 @@ enum UpdateFlags { * Update the second derivatives of the * shape functions on the real cell. */ - update_second_derivatives + update_second_derivatives = 64 }; diff --git a/deal.II/deal.II/include/fe/fe_values.h b/deal.II/deal.II/include/fe/fe_values.h index 303c126a24..35c15b2f76 100644 --- a/deal.II/deal.II/include/fe/fe_values.h +++ b/deal.II/deal.II/include/fe/fe_values.h @@ -148,9 +148,9 @@ template class Quadrature; * a reference to a whole field. Usually these fields contain * the values of all trial functions at all quadrature points. * - * \item #get_function_values#, #get_function_gradients#: these - * two functions offer a simple way to avoid the detour of the - * trial functions, if you have a finite solution (resp. the + * \item #get_function_values#, #get_function_grads#, #...#: these + * functions offer a simple way to avoid the detour of the + * trial functions, if you have a finite element solution (resp. the * vector of values associated with the different trial functions.) * Then you may want to get information from the restriction of * the finite element function to a certain cell, e.g. the values @@ -159,7 +159,8 @@ template class Quadrature; * you pass it a vector holding the finite element solution and the * functions return the values or gradients of the finite element * function restricted to the cell which was given last time the - * #reinit# function was given. + * #reinit# function was given. The same applies for the functions + * returning higher derivatives of the solution. * * Though possible in principle, these functions do not call the * #reinit# function, you have to do so yourself beforehand. On the @@ -316,7 +317,7 @@ class FEValuesBase { * documentation for the matrix itself. */ const vector > > & get_shape_grads () const; - + /** * Return the gradients of the finite * element function characterized @@ -330,6 +331,53 @@ class FEValuesBase { void get_function_grads (const dVector &fe_function, vector > &gradients) const; + /** + * Return the 2nd derivatives of + * the #i#th shape function at + * the #j# quadrature point. 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 derivatives' + * values is returned, there + * should be no major performance + * drawback. The function + * returns the derivatives on the + * real element, not the + * reference element. + */ + const Tensor<2,dim> & shape_2nd_derivative (const unsigned int function, + const unsigned int quadrature_point) const; + + /** + * Return a pointer to the + * matrix holding all 2nd + * derivatives of shape functions + * at all integration points, on + * the present cell. For the + * format of this matrix, see the + * documentation for the matrix + * itself. + */ + const vector > > & get_shape_2nd_derivatives () const; + + /** + * Return the tensor of second + * derivatives of the finite + * element function characterized + * by #fe_function# restricted to + * #cell# at the quadrature points. + * + * The function assumes that the + * #second_derivatives# object already has + * the right size. + */ + void get_function_2nd_derivatives (const dVector &fe_function, + vector > &second_derivatives) const; + /** * Return the position of the #i#th * quadrature point in real space. @@ -1155,6 +1203,16 @@ FEValuesBase::get_shape_grads () const { +template +inline +const vector > > & +FEValuesBase::get_shape_2nd_derivatives () const { + Assert (update_flags & update_second_derivatives, ExcAccessToUninitializedField()); + return shape_2nd_derivatives; +}; + + + template inline const vector > & diff --git a/deal.II/deal.II/source/fe/fe.cc b/deal.II/deal.II/source/fe/fe.cc index 57151da0e9..7fe883f092 100644 --- a/deal.II/deal.II/source/fe/fe.cc +++ b/deal.II/deal.II/source/fe/fe.cc @@ -318,6 +318,8 @@ void FiniteElement::fill_fe_values (const DoFHandler::cell_iterator &, const vector > &, vector > &, const bool, + vector > &, + const bool, vector > &, const bool, vector > &, @@ -337,6 +339,8 @@ void FiniteElement::fill_fe_face_values (const DoFHandler::cell_iterat const vector > &global_unit_points, vector > &jacobians, const bool compute_jacobians, + vector > &jacobians_grad, + const bool compute_jacobians_grad, vector > &support_points, const bool compute_support_points, vector > &q_points, @@ -360,6 +364,7 @@ void FiniteElement::fill_fe_face_values (const DoFHandler::cell_iterat vector > dummy(total_dofs); fill_fe_values (cell, global_unit_points, jacobians, compute_jacobians, + jacobians_grad, compute_jacobians_grad, dummy, false, q_points, compute_q_points, shape_values_transform, shape_gradients_transform, @@ -388,6 +393,8 @@ void FiniteElement::fill_fe_subface_values (const DoFHandler::cell_ite const vector > &global_unit_points, vector > &jacobians, const bool compute_jacobians, + vector > &jacobians_grad, + const bool compute_jacobians_grad, vector > &q_points, const bool compute_q_points, vector &face_jacobi_determinants, @@ -407,6 +414,7 @@ void FiniteElement::fill_fe_subface_values (const DoFHandler::cell_ite vector > dummy(total_dofs); fill_fe_values (cell, global_unit_points, jacobians, compute_jacobians, + jacobians_grad, compute_jacobians_grad, dummy, false, q_points, compute_q_points, shape_values_transform, shape_gradients_transform, diff --git a/deal.II/deal.II/source/fe/fe_lib.criss_cross.cc b/deal.II/deal.II/source/fe/fe_lib.criss_cross.cc index 0af836516f..de29c860ed 100644 --- a/deal.II/deal.II/source/fe/fe_lib.criss_cross.cc +++ b/deal.II/deal.II/source/fe/fe_lib.criss_cross.cc @@ -928,12 +928,14 @@ template void FECrissCross::fill_fe_values (const DoFHandler::cell_iterator &cell, const vector > &unit_points, vector > &jacobians, - const bool compute_jacobians, - vector > &support_points, + const bool compute_jacobians, + vector > &jacobians_grad, + const bool compute_jacobians_grad, + vector > &support_points, const bool, - vector > &q_points, - const bool compute_q_points, - const dFMatrix &shape_values_transform, + vector > &q_points, + const bool compute_q_points, + const dFMatrix &shape_values_transform, const vector > > &/*shape_grad_transform*/, const Boundary &boundary) const { Assert (jacobians.size() == unit_points.size(), @@ -1051,6 +1053,92 @@ void FECrissCross::fill_fe_values (const DoFHandler::cell_iterator &ce // not implemented at present Assert (false, ExcNotImplemented()); }; + + + if (compute_jacobians_grad) + switch (dim) + { + case 1: + { + // derivative of the + // jacobian is always zero + // for a linear mapping in 1d + for (unsigned int point=0; point void FELinearMapping::fill_fe_values (const DoFHandler::cell_iterator &cell, const vector > &unit_points, vector > &jacobians, - const bool compute_jacobians, + const bool compute_jacobians, + vector > &jacobians_grad, + const bool compute_jacobians_grad, vector > &support_points, const bool compute_support_points, vector > &q_points, @@ -379,7 +381,8 @@ void FELinearMapping::fill_fe_values (const DoFHandler::cell_iterator is multiplied to the unit cell gradient *from the right*! be very careful with these things. - The following little program tests the correct behaviour: + The following little program tests the correct behaviour. The program can + also be found in the directory. ------------------------------------------- #include @@ -488,7 +491,94 @@ void FELinearMapping::fill_fe_values (const DoFHandler::cell_iterator // not implemented at present Assert (false, ExcNotImplemented()); }; + + if (compute_jacobians_grad) + switch (dim) + { + case 1: + { + // derivative of the + // jacobian is always zero + // for a linear mapping in 1d + for (unsigned int point=0; point::FEValuesBase (const unsigned int n_q_points, n_transform_functions (n_transform_functions), shape_values (n_values_arrays, dFMatrix(n_dofs, n_q_points)), shape_gradients (n_dofs, vector >(n_q_points)), + shape_2nd_derivatives (n_dofs, vector >(n_q_points)), weights (n_q_points, 0), JxW_values (n_q_points, 0), quadrature_points (n_q_points, Point()), @@ -127,6 +128,48 @@ void FEValuesBase::get_function_grads (const dVector &fe_function, +template +const Tensor<2,dim> & +FEValuesBase::shape_2nd_derivative (const unsigned int i, + const unsigned int j) const { + Assert (i +void FEValuesBase::get_function_2nd_derivatives (const dVector &fe_function, + vector > &second_derivatives) const { + Assert (second_derivatives.size() == n_quadrature_points, + ExcWrongVectorSize(second_derivatives.size(), n_quadrature_points)); + + // get function values of dofs + // on this cell + dVector dof_values (total_dofs); + present_cell->get_dof_values (fe_function, dof_values); + + // initialize with zero + fill_n (second_derivatives.begin(), n_quadrature_points, Tensor<2,dim>()); + + // add up contributions of trial + // functions + for (unsigned int point=0; point tmp(shape_2nd_derivatives[shape_func][point]); + tmp *= dof_values(shape_func); + second_derivatives[point] += tmp; + }; +}; + + + template const Point & FEValuesBase::quadrature_point (const unsigned int i) const { @@ -218,10 +261,11 @@ void FEValues::reinit (const typename DoFHandler::cell_iterator &cell, present_cell = cell; // fill jacobi matrices and real // quadrature points - if ((update_flags & update_jacobians) || - (update_flags & update_JxW_values)|| - (update_flags & update_q_points) || - (update_flags & update_gradients) || + if ((update_flags & update_jacobians) || + (update_flags & update_JxW_values) || + (update_flags & update_q_points) || + (update_flags & update_gradients) || + (update_flags & update_second_derivatives) || (update_flags & update_support_points)) fe->fill_fe_values (cell, unit_quadrature_points, @@ -230,8 +274,8 @@ void FEValues::reinit (const typename DoFHandler::cell_iterator &cell, update_JxW_values | update_gradients | update_second_derivatives), -// jacobi_matrices_grad, -// update_flags & update_second_derivatives, + jacobi_matrices_grad, + update_flags & update_second_derivatives, support_points, update_flags & update_support_points, quadrature_points, @@ -407,31 +451,35 @@ void FEFaceValues::reinit (const typename DoFHandler::cell_iterator &c selected_dataset = face_no; // fill jacobi matrices and real // quadrature points - if ((update_flags & update_jacobians) || - (update_flags & update_JxW_values)|| - (update_flags & update_q_points) || - (update_flags & update_gradients) || - (update_flags & update_support_points) || + if ((update_flags & update_jacobians) || + (update_flags & update_JxW_values) || + (update_flags & update_q_points) || + (update_flags & update_gradients) || + (update_flags & update_second_derivatives) || + (update_flags & update_support_points) || (update_flags & update_JxW_values)) fe->fill_fe_face_values (cell, - face_no, - unit_face_quadrature_points, - unit_quadrature_points[face_no], - jacobi_matrices, - update_flags & (update_jacobians | - update_gradients | - update_JxW_values), - support_points, - update_flags & update_support_points, - quadrature_points, - update_flags & update_q_points, - face_jacobi_determinants, - update_flags & update_JxW_values, - normal_vectors, - update_flags & update_normal_vectors, - shape_values_transform[face_no], - unit_shape_gradients_transform[face_no], - boundary); + face_no, + unit_face_quadrature_points, + unit_quadrature_points[face_no], + jacobi_matrices, + update_flags & (update_jacobians | + update_gradients | + update_JxW_values | + update_second_derivatives), + jacobi_matrices_grad, + update_flags & update_second_derivatives, + support_points, + update_flags & update_support_points, + quadrature_points, + update_flags & update_q_points, + face_jacobi_determinants, + update_flags & update_JxW_values, + normal_vectors, + update_flags & update_normal_vectors, + shape_values_transform[face_no], + unit_shape_gradients_transform[face_no], + boundary); // compute gradients on real element if // requested @@ -579,29 +627,33 @@ void FESubfaceValues::reinit (const typename DoFHandler::cell_iterator selected_dataset = face_no*(1<<(dim-1)) + subface_no; // fill jacobi matrices and real // quadrature points - if ((update_flags & update_jacobians) || - (update_flags & update_JxW_values)|| - (update_flags & update_q_points) || - (update_flags & update_gradients) || + if ((update_flags & update_jacobians) || + (update_flags & update_JxW_values) || + (update_flags & update_q_points) || + (update_flags & update_gradients) || + (update_flags & update_second_derivatives) || (update_flags & update_JxW_values)) fe->fill_fe_subface_values (cell, - face_no, - subface_no, - unit_face_quadrature_points, - unit_quadrature_points[selected_dataset], - jacobi_matrices, - update_flags & (update_jacobians | - update_gradients | - update_JxW_values), - quadrature_points, - update_flags & update_q_points, - face_jacobi_determinants, - update_flags & update_JxW_values, - normal_vectors, - update_flags & update_normal_vectors, - shape_values_transform[selected_dataset], - unit_shape_gradients_transform[selected_dataset], - boundary); + face_no, + subface_no, + unit_face_quadrature_points, + unit_quadrature_points[selected_dataset], + jacobi_matrices, + update_flags & (update_jacobians | + update_gradients | + update_JxW_values| + update_second_derivatives), + jacobi_matrices_grad, + update_flags & update_second_derivatives, + quadrature_points, + update_flags & update_q_points, + face_jacobi_determinants, + update_flags & update_JxW_values, + normal_vectors, + update_flags & update_normal_vectors, + shape_values_transform[selected_dataset], + unit_shape_gradients_transform[selected_dataset], + boundary); // compute gradients on real element if // requested diff --git a/deal.II/deal.II/source/fe/scripts/2d/postprocess b/deal.II/deal.II/source/fe/scripts/2d/postprocess index 57425e87fc..e7034ba430 100644 --- a/deal.II/deal.II/source/fe/scripts/2d/postprocess +++ b/deal.II/deal.II/source/fe/scripts/2d/postprocess @@ -22,6 +22,7 @@ perl -pi -e 's/\[(\d)\]\[(\d)\] =/($1,$2) =/g;' *2d_inverse_jacobian perl -pi -e 's/^\s*t/const double t/g;' *2d_inverse_jacobian_grad perl -pi -e 's/x\[(\d)\]/vertices[$1](0)/g;' *2d_inverse_jacobian_grad perl -pi -e 's/y\[(\d)\]/vertices[$1](1)/g;' *2d_inverse_jacobian_grad +perl -pi -e 's/inverseJacobian/jacobians_grad[point]/g;' *2d_inverse_jacobian_grad perl -pi -e 's/([^;])\n/$1/g;' *2d_shape_grad perl -pi -e 's/grad_phi_polynom\[(\d+)\]\[0\] = (.*);/case $1: return Point<2>($2,/g;' *2d_shape_grad