From: Martin Kronbichler Date: Sat, 8 Aug 2015 20:28:39 +0000 (+0200) Subject: Also compute gradients of Jacobians on faces X-Git-Tag: v8.4.0-rc2~641^2 X-Git-Url: https://gitweb.dealii.org/cgi-bin/gitweb.cgi?a=commitdiff_plain;h=refs%2Fpull%2F1311%2Fhead;p=dealii.git Also compute gradients of Jacobians on faces --- diff --git a/doc/news/changes.h b/doc/news/changes.h index b687662028..7491822f11 100644 --- a/doc/news/changes.h +++ b/doc/news/changes.h @@ -90,6 +90,13 @@ inconvenience this causes.
    +
  1. New: FEFaceValues and FESubfaceValues can now also compute + gradients of the Jacobian of the transformation from unit to real cell, + controlled by update_jacobian_grads. +
    + (Martin Kronbichler, 2015/08/08) +
    2. +
  2. New: There is now a function MemoryConsumption::memory_consumption() for std_cxx11::unique_ptr arguments.
    diff --git a/source/fe/fe_values.cc b/source/fe/fe_values.cc index 6f18dc0fbe..a5a19d25fb 100644 --- a/source/fe/fe_values.cc +++ b/source/fe/fe_values.cc @@ -3709,9 +3709,6 @@ void FEFaceValues::do_reinit (const unsigned int face_no) const typename Triangulation::cell_iterator cell=*this->present_cell; this->present_face_index=cell->face_index(face_no); - Assert(!(this->update_flags & update_jacobian_grads), - ExcNotImplemented()); - this->get_mapping().fill_fe_face_values(*this->present_cell, face_no, this->quadrature, @@ -3943,9 +3940,6 @@ void FESubfaceValues::do_reinit (const unsigned int face_no, this->present_face_index=subface_index; } - Assert(!(this->update_flags & update_jacobian_grads), - ExcNotImplemented()); - // now ask the mapping and the finite element to do the actual work this->get_mapping().fill_fe_subface_values(*this->present_cell, face_no, diff --git a/source/fe/mapping_q1.cc b/source/fe/mapping_q1.cc index 4d9bde09b6..58cfd8ceb3 100644 --- a/source/fe/mapping_q1.cc +++ b/source/fe/mapping_q1.cc @@ -835,6 +835,54 @@ namespace internal } } + + /** + * Update the Hessian of the transformation from unit to real cell, the + * Jacobian gradients. + * + * Skip the computation if possible as indicated by the first argument. + */ + template + void + maybe_update_jacobian_grads (const CellSimilarity::Similarity cell_similarity, + const typename dealii::MappingQ1::DataSetDescriptor data_set, + const typename dealii::MappingQ1::InternalData &data, + std::vector > &jacobian_grads) + { + const UpdateFlags update_flags(data.current_update_flags()); + if (update_flags & update_jacobian_grads) + { + const unsigned int n_q_points = jacobian_grads.size(); + + if (cell_similarity != CellSimilarity::translation) + { + for (unsigned int point=0; point *second = + &data.second_derivative(point+data_set, 0); + double result [spacedim][dim][dim]; + for (unsigned int i=0; i::cell_iterator &cell, output_data.jacobians[point] = data.contravariant[point]; } - // calculate values of the derivatives of the Jacobians. do it here, since - // we only do it for cells, not faces. - if (update_flags & update_jacobian_grads) - { - AssertDimension (output_data.jacobian_grads.size(), n_q_points); - - if (cell_similarity != CellSimilarity::translation) - { - - std::fill(output_data.jacobian_grads.begin(), - output_data.jacobian_grads.end(), - DerivativeForm<2,dim,spacedim>()); - - - const unsigned int data_set = DataSetDescriptor::cell(); - - for (unsigned int point=0; point *second = - &data.second_derivative(point+data_set, 0); - double result [spacedim][dim][dim]; - for (unsigned int i=0; i::cell_iterator &cell, output_data.inverse_jacobians[point] = data.covariant[point].transpose(); } + internal::maybe_update_jacobian_grads (cell_similarity, + DataSetDescriptor::cell (), + data, + output_data.jacobian_grads); + return cell_similarity; } @@ -1220,6 +1229,10 @@ namespace internal maybe_update_Jacobians (CellSimilarity::none, data_set, data); + maybe_update_jacobian_grads (CellSimilarity::none, + data_set, + data, + output_data.jacobian_grads); maybe_compute_face_data (mapping, cell, face_no, subface_no, quadrature.size(), quadrature.get_weights(), data, diff --git a/tests/fe/jacobians.cc b/tests/fe/jacobians.cc index f9427a0c03..696de23412 100644 --- a/tests/fe/jacobians.cc +++ b/tests/fe/jacobians.cc @@ -14,8 +14,8 @@ // --------------------------------------------------------------------- -// Show the Jacobians and inverse Jacobians on hyperball with one quadrature -// point +// Show the Jacobians, inverse Jacobians, and Jacobian gradients on hyperball +// with one quadrature point #include "../tests.h" #include @@ -44,7 +44,8 @@ void test() quad_p(d) = 0.42 + 0.11 * d; Quadrature quad(quad_p); FEValues fe_val (mapping, dummy, quad, - update_jacobians | update_inverse_jacobians); + update_jacobians | update_inverse_jacobians | + update_jacobian_grads); deallog << dim << "d Jacobians:" << std::endl; typename Triangulation::active_cell_iterator cell = tria.begin_active(), endc = tria.end(); @@ -72,6 +73,21 @@ void test() deallog << std::endl; } deallog << std::endl; + + deallog << dim << "d Jacobian gradients:" << std::endl; + cell = tria.begin_active(); + endc = tria.end(); + for ( ; cell != endc; ++cell) + { + fe_val.reinit (cell); + + for (unsigned int d=0; d quad(quad_p); FEFaceValues fe_val (mapping, dummy, quad, - update_jacobians | update_inverse_jacobians); + update_jacobians | update_inverse_jacobians | + update_jacobian_grads); FESubfaceValues fe_sub_val (mapping, dummy, quad, - update_jacobians | update_inverse_jacobians); + update_jacobians | update_inverse_jacobians | + update_jacobian_grads); deallog << dim << "d Jacobians:" << std::endl; typename Triangulation::active_cell_iterator cell = tria.begin_active(), endc = tria.end(); @@ -104,6 +106,36 @@ void test() } deallog << std::endl; + + deallog << dim << "d Jacobian gradients:" << std::endl; + cell = tria.begin_active(); + endc = tria.end(); + for ( ; cell != endc; ++cell) + for (unsigned int f=0; f::faces_per_cell; ++f) + { + fe_val.reinit (cell, f); + + for (unsigned int d=0; dat_boundary(f) == false && + cell->neighbor(f)->level() < cell->level()) + { + fe_sub_val.reinit(cell->neighbor(f), cell->neighbor_face_no(f), + cell->neighbor_of_coarser_neighbor(f).second); + + for (unsigned int d=0; d @@ -44,9 +44,11 @@ void test() quad_p(d) = 0.42 + 0.11 * d; Quadrature quad(quad_p); FEFaceValues fe_val (mapping, dummy, quad, - update_jacobians | update_inverse_jacobians); + update_jacobians | update_inverse_jacobians | + update_jacobian_grads); FESubfaceValues fe_sub_val (mapping, dummy, quad, - update_jacobians | update_inverse_jacobians); + update_jacobians | update_inverse_jacobians | + update_jacobian_grads); deallog << dim << "d Jacobians:" << std::endl; typename Triangulation::active_cell_iterator cell = tria.begin_active(), endc = tria.end(); @@ -103,6 +105,36 @@ void test() } deallog << std::endl; + + deallog << dim << "d Jacobian gradients:" << std::endl; + cell = tria.begin_active(); + endc = tria.end(); + for ( ; cell != endc; ++cell) + for (unsigned int f=0; f::faces_per_cell; ++f) + { + fe_val.reinit (cell, f); + + for (unsigned int d=0; dat_boundary(f) == false && + cell->neighbor(f)->level() < cell->level()) + { + fe_sub_val.reinit(cell->neighbor(f), cell->neighbor_face_no(f), + cell->neighbor_of_coarser_neighbor(f).second); + + for (unsigned int d=0; d