From 52c4c7dd171110845a7e682318ff0e6218d5326d Mon Sep 17 00:00:00 2001 From: bangerth Date: Thu, 13 May 2010 13:22:15 +0000 Subject: [PATCH] Patch from Markus Buerg: Implement FEValuesViews::Vector::curl and friends. git-svn-id: https://svn.dealii.org/trunk@21127 0785d39b-7218-0410-832d-ea1e28bc413d --- deal.II/deal.II/include/fe/fe_values.h | 163 ++++++++++++++++++++++++- deal.II/deal.II/source/fe/fe_values.cc | 161 +++++++++++++++++++++++- deal.II/doc/news/changes.h | 7 ++ 3 files changed, 329 insertions(+), 2 deletions(-) diff --git a/deal.II/deal.II/include/fe/fe_values.h b/deal.II/deal.II/include/fe/fe_values.h index 5d1b482199..5eec718392 100644 --- a/deal.II/deal.II/include/fe/fe_values.h +++ b/deal.II/deal.II/include/fe/fe_values.h @@ -500,6 +500,19 @@ namespace FEValuesViews */ typedef double divergence_type; + /** + * A typedef for the type of the curl + * of the view this class + * represents. Here, for a set of + * dim components of the + * finite element, the curl is of type + * Tensor@<1, 1@> (i.e. a + * scalar) in 2d, and of type + * Tensor@<1, 3@> (i.e. a + * vector) in 3d. + */ + typedef Tensor<1, (spacedim == 3)? 3 : 1> curl_type; + /** * A typedef for the type of second * derivatives of the view this class @@ -598,6 +611,17 @@ namespace FEValuesViews divergence (const unsigned int shape_function, const unsigned int q_point) const; + /** + * Return the vector curl of + * the vector components selected by + * this view, for the shape function + * and quadrature point selected by the + * arguments. + */ + curl_type + curl (const unsigned int shape_function, + const unsigned int q_point) const; + /** * Return the Hessian (the tensor of * rank 2 of all second derivatives) of @@ -689,7 +713,7 @@ namespace FEValuesViews * * There is no equivalent function such * as - * FEValuesBase::get_function_gradients + * FEValuesBase::get_function_divergences * in the FEValues classes but the * information can be obtained from * FEValuesBase::get_function_gradients, @@ -699,6 +723,29 @@ namespace FEValuesViews void get_function_divergences (const InputVector& fe_function, std::vector& divergences) const; + /** + * Return the curl of the selected + * vector components of the finite + * element function characterized by + * fe_function at the + * quadrature points of the cell, face + * or subface selected the last time + * the reinit function of the + * FEValues object was called, at the + * quadrature points. + * + * There is no equivalent function such + * as + * FEValuesBase::get_function_curls + * in the FEValues classes but the + * information can be obtained from + * FEValuesBase::get_function_gradients, + * of course. + */ + template + void get_function_curls (const InputVector& fe_function, + std::vector& curls) const; + /** * Return the Hessians of the selected * vector components of the finite @@ -3730,6 +3777,120 @@ namespace FEValuesViews template inline + typename Vector::curl_type + Vector::curl (const unsigned int shape_function, + const unsigned int q_point) const + { + // this function works like in the case + // above + typedef FEValuesBase FVB; + + Assert (shape_function < fe_values.fe->dofs_per_cell, + ExcIndexRange (shape_function, 0, fe_values.fe->dofs_per_cell)); + Assert (fe_values.update_flags & update_gradients, + typename FVB::ExcAccessToUninitializedField()); + // same as for the scalar case except + // that we have one more index + const int snc = shape_function_data[shape_function].single_nonzero_component; + + if (snc == -2) + return curl_type (); + + else + switch (dim) { + case 1: + return curl_type (); + + case 2: { + if (snc != -1) { + curl_type return_value; + + switch (shape_function_data[shape_function].single_nonzero_component_index) { + case 0: { + return_value[0] = -1.0 * fe_values.shape_gradients[snc][q_point][1]; + return return_value; + } + + default: { + return_value[0] = fe_values.shape_gradients[snc][q_point][2]; + return return_value; + } + } + } + + else { + curl_type return_value; + + return_value[0] = 0.0; + + if (shape_function_data[shape_function].is_nonzero_shape_function_component[0]) + return_value[0] -= fe_values.shape_gradients[shape_function_data[shape_function].row_index[0]][q_point][1]; + + if (shape_function_data[shape_function].is_nonzero_shape_function_component[1]) + return_value[0] += fe_values.shape_gradients[shape_function_data[shape_function].row_index[1]][q_point][0]; + + return return_value; + } + } + + case 3: { + if (snc != -1) { + curl_type return_value; + + switch (shape_function_data[shape_function].single_nonzero_component_index) { + case 0: { + return_value[0] = 0; + return_value[1] = fe_values.shape_gradients[snc][q_point][2]; + return_value[2] = -1.0 * fe_values.shape_gradients[snc][q_point][1]; + return return_value; + } + + case 1: { + return_value[0] = -1.0 * fe_values.shape_gradients[snc][q_point][2]; + return_value[1] = 0; + return_value[2] = fe_values.shape_gradients[snc][q_point][0]; + return return_value; + } + + default: { + return_value[0] = fe_values.shape_gradients[snc][q_point][1]; + return_value[1] = -1.0 * fe_values.shape_gradients[snc][q_point][0]; + return_value[2] = 0; + return return_value; + } + } + } + + else { + curl_type return_value; + + for (unsigned int i = 0; i < dim; ++i) + return_value[i] = 0.0; + + if (shape_function_data[shape_function].is_nonzero_shape_function_component[0]) { + return_value[1] += fe_values.shape_gradients[shape_function_data[shape_function].row_index[0]][q_point][2]; + return_value[2] -= fe_values.shape_gradients[shape_function_data[shape_function].row_index[0]][q_point][1]; + } + + if (shape_function_data[shape_function].is_nonzero_shape_function_component[1]) { + return_value[0] -= fe_values.shape_gradients[shape_function_data[shape_function].row_index[1]][q_point][2]; + return_value[2] += fe_values.shape_gradients[shape_function_data[shape_function].row_index[1]][q_point][0]; + } + + if (shape_function_data[shape_function].is_nonzero_shape_function_component[2]) { + return_value[0] += fe_values.shape_gradients[shape_function_data[shape_function].row_index[2]][q_point][1]; + return_value[1] -= fe_values.shape_gradients[shape_function_data[shape_function].row_index[2]][q_point][0]; + } + + return return_value; + } + } + } + } + + + template + inline typename Vector::hessian_type Vector::hessian (const unsigned int shape_function, const unsigned int q_point) const diff --git a/deal.II/deal.II/source/fe/fe_values.cc b/deal.II/deal.II/source/fe/fe_values.cc index c141fe2771..ff01035778 100644 --- a/deal.II/deal.II/source/fe/fe_values.cc +++ b/deal.II/deal.II/source/fe/fe_values.cc @@ -755,9 +755,168 @@ namespace FEValuesViews } } + template + template + void + Vector::get_function_curls (const InputVector &fe_function, + std::vector &curls) const { + typedef FEValuesBase FVB; + + Assert (fe_values.update_flags & update_gradients, + typename FVB::ExcAccessToUninitializedField()); + Assert (curls.size() == fe_values.n_quadrature_points, + ExcDimensionMismatch (curls.size(), fe_values.n_quadrature_points)); + Assert (fe_values.present_cell.get () != 0, + ExcMessage ("FEValues object is not reinit'ed to any cell")); + Assert (fe_function.size () == fe_values.present_cell->n_dofs_for_dof_handler (), + ExcDimensionMismatch (fe_function.size (), fe_values.present_cell->n_dofs_for_dof_handler ())); + // get function values of dofs on this + // cell + dealii::Vector dof_values (fe_values.dofs_per_cell); + fe_values.present_cell->get_interpolated_dof_values (fe_function, dof_values); + + std::fill (curls.begin (), curls.end (), curl_type ()); + + switch (dim) { + case 2: { + for (unsigned int shape_function = 0; shape_function < fe_values.fe->dofs_per_cell; ++shape_function) { + const int snc = shape_function_data[shape_function].single_nonzero_component; + + if (snc == -2) + // shape function is zero for the selected components + continue; + + const double value = dof_values (shape_function); + + if (value == 0.) + continue; + + if (snc != -1) { + const unsigned int comp = shape_function_data[shape_function].single_nonzero_component_index; + const Tensor<1, spacedim> *shape_gradient_ptr = &fe_values.shape_gradients[snc][0]; + + switch (shape_function_data[shape_function].single_nonzero_component_index) { + case 0: { + for (unsigned int q_point = 0; q_point < fe_values.n_quadrature_points; ++q_point) + curls[q_point][0] = -1.0 * value * (*shape_gradient_ptr++)[1]; + + break; + } + + default: + for (unsigned int q_point = 0; q_point < fe_values.n_quadrature_points; ++q_point) + curls[q_point][0] = value * (*shape_gradient_ptr)[0]; + } + } + + else { + for (unsigned int q_point = 0; q_point < fe_values.n_quadrature_points; ++q_point) + curls[q_point][0] = 0; + + if (shape_function_data[shape_function].is_nonzero_shape_function_component[0]) { + const Tensor<1,spacedim> *shape_gradient_ptr = &fe_values.shape_gradients[shape_function_data[shape_function].row_index[0]][0]; + + for (unsigned int q_point = 0; q_point < fe_values.n_quadrature_points; ++q_point) + curls[q_point][0] -= value * (*shape_gradient_ptr++)[1]; + } + + if (shape_function_data[shape_function].is_nonzero_shape_function_component[1]) { + const Tensor<1,spacedim> *shape_gradient_ptr = &fe_values.shape_gradients[shape_function_data[shape_function].row_index[1]][0]; + + for (unsigned int q_point = 0; q_point < fe_values.n_quadrature_points; ++q_point) + curls[q_point][0] += value * (*shape_gradient_ptr++)[0]; + } + } + } + break; + } + + case 3: { + for (unsigned int shape_function = 0; shape_function < fe_values.fe->dofs_per_cell; ++shape_function) { + const int snc = shape_function_data[shape_function].single_nonzero_component; + + if (snc == -2) + // shape function is zero for the selected components + continue; + + const double value = dof_values (shape_function); + + if (value == 0.) + continue; + + if (snc != -1) { + const unsigned int comp = shape_function_data[shape_function].single_nonzero_component_index; + const Tensor<1, spacedim> *shape_gradient_ptr = &fe_values.shape_gradients[snc][0]; + + switch (shape_function_data[shape_function].single_nonzero_component_index) { + case 0: { + for (unsigned int q_point = 0; q_point < fe_values.n_quadrature_points; ++q_point) { + curls[q_point][0] = 0; + curls[q_point][1] = value * (*shape_gradient_ptr)[2]; + curls[q_point][2] = -1.0 * value * (*shape_gradient_ptr++)[1]; + } + + break; + } + + case 1: { + for (unsigned int q_point = 0; q_point < fe_values.n_quadrature_points; ++q_point) { + curls[q_point][0] = -1.0 * value * (*shape_gradient_ptr)[2]; + curls[q_point][1] = 0; + curls[q_point][2] = value * (*shape_gradient_ptr++)[0]; + } + + break; + } + + default: + for (unsigned int q_point = 0; q_point < fe_values.n_quadrature_points; ++q_point) { + curls[q_point][0] = value * (*shape_gradient_ptr)[1]; + curls[q_point][1] = -1.0 * value * (*shape_gradient_ptr++)[0]; + curls[q_point][2] = 0; + } + } + } + + else { + for (unsigned int q_point = 0; q_point < fe_values.n_quadrature_points; ++q_point) + for (unsigned int d = 0; d < dim; ++d) + curls[q_point][d] = 0; + + if (shape_function_data[shape_function].is_nonzero_shape_function_component[0]) { + const Tensor<1,spacedim> *shape_gradient_ptr = &fe_values.shape_gradients[shape_function_data[shape_function].row_index[0]][0]; + + for (unsigned int q_point = 0; q_point < fe_values.n_quadrature_points; ++q_point) { + curls[q_point][1] += value * (*shape_gradient_ptr++)[2]; + curls[q_point][2] -= value * (*shape_gradient_ptr++)[1]; + } + } + + if (shape_function_data[shape_function].is_nonzero_shape_function_component[1]) { + const Tensor<1,spacedim> *shape_gradient_ptr = &fe_values.shape_gradients[shape_function_data[shape_function].row_index[1]][0]; + + for (unsigned int q_point = 0; q_point < fe_values.n_quadrature_points; ++q_point) { + curls[q_point][0] -= value * (*shape_gradient_ptr++)[2]; + curls[q_point][2] += value * (*shape_gradient_ptr++)[0]; + } + } + + if (shape_function_data[shape_function].is_nonzero_shape_function_component[2]) { + const Tensor<1,spacedim> *shape_gradient_ptr = &fe_values.shape_gradients[shape_function_data[shape_function].row_index[2]][0]; + + for (unsigned int q_point = 0; q_point < fe_values.n_quadrature_points; ++q_point) { + curls[q_point][0] += value * (*shape_gradient_ptr++)[1]; + curls[q_point][1] -= value * (*shape_gradient_ptr++)[0]; + } + } + } + } + } + } + } - template + template template void Vector:: diff --git a/deal.II/doc/news/changes.h b/deal.II/doc/news/changes.h index cae2c5aa8c..e216d88a25 100644 --- a/deal.II/doc/news/changes.h +++ b/deal.II/doc/news/changes.h @@ -636,6 +636,13 @@ inconvenience this causes.

deal.II

    +
  1. +

    New: The FEValuesViews::Vector class now has functions + FEValuesViews::Vector::curl and FEValuesViews::Vector::get_function_curls. +
    + (Markus Bürg 2010/05/13) +

    2. +

  2. New: TriaAccessor::extent_in_direction() returns the length of an object in a given direction. -- 2.39.5