From 483d5cdc3b9bf791736945fa4bc194f78fc29528 Mon Sep 17 00:00:00 2001 From: leicht Date: Wed, 17 Jan 2007 15:38:04 +0000 Subject: [PATCH] reorder shape functions on faces of 3D cells if the face has non-standard face_orientation, currently only works for FE_Q git-svn-id: https://svn.dealii.org/trunk@14330 0785d39b-7218-0410-832d-ea1e28bc413d --- deal.II/deal.II/include/fe/fe.h | 25 ++++ deal.II/deal.II/include/fe/fe_q.h | 20 ++++ deal.II/deal.II/include/fe/fe_system.h | 19 +++ deal.II/deal.II/include/fe/fe_values.h | 157 +++++++++++++++++++------ deal.II/deal.II/source/fe/fe.cc | 24 ++++ deal.II/deal.II/source/fe/fe_q.cc | 31 +++++ deal.II/deal.II/source/fe/fe_system.cc | 34 ++++++ deal.II/deal.II/source/fe/fe_values.cc | 72 ++++++++++++ deal.II/doc/news/changes.html | 10 ++ 9 files changed, 356 insertions(+), 36 deletions(-) diff --git a/deal.II/deal.II/include/fe/fe.h b/deal.II/deal.II/include/fe/fe.h index 6b1eb6e037..60107d9021 100644 --- a/deal.II/deal.II/include/fe/fe.h +++ b/deal.II/deal.II/include/fe/fe.h @@ -1459,6 +1459,31 @@ class FiniteElement : public Subscriptor, unsigned int component_to_block_index (const unsigned int component) const; + /** + * For faces with non-standard + * face_orientation in 3D, the shape + * functions on faces have to be permutated + * in order to be combined with the correct + * dofs. This function returns a vector of + * integer values, that have to be added to + * the index of a shape function in order + * to result in the permutated index. Prior + * content of the vector @p shifts is + * erased. In 3D a vector of length @p + * dofs_per_quad is returned, in 2D and 1D + * there is no need for permutation and a + * vector of length 0 is returned, the same + * is true for elements which have no dofs + * on quads. The general implementation + * returns a vector of zeros, resulting in + * no permutation at all. This has to be + * overloaded by derived finite element + * classes. + */ + virtual + void + get_face_shape_function_shifts (std::vector &shifts) const; + //@} /** diff --git a/deal.II/deal.II/include/fe/fe_q.h b/deal.II/deal.II/include/fe/fe_q.h index 5babd3ffed..92d5ab0a1e 100644 --- a/deal.II/deal.II/include/fe/fe_q.h +++ b/deal.II/deal.II/include/fe/fe_q.h @@ -438,6 +438,26 @@ class FE_Q : public FE_Poly,dim> virtual FiniteElementDomination::Domination compare_for_face_domination (const FiniteElement &fe_other) const; + + /** + * For faces with non-standard + * face_orientation in 3D, the shape + * functions on faces have to be permutated + * in order to be combined with the correct + * dofs. This function returns a vector of + * integer values, that have to be added to + * the index of a shape function in order + * to result in the permutated index. Prior + * content of the vector @p shifts is + * erased. In 3D a vector of length @p + * dofs_per_quad is returned, in 2D and 1D + * there is no need for permutation and a + * vector of length 0 is returned. + */ + virtual + void + get_face_shape_function_shifts (std::vector &shifts) const; + //@} /** diff --git a/deal.II/deal.II/include/fe/fe_system.h b/deal.II/deal.II/include/fe/fe_system.h index 8968067916..d2d8ac1c65 100644 --- a/deal.II/deal.II/include/fe/fe_system.h +++ b/deal.II/deal.II/include/fe/fe_system.h @@ -539,6 +539,25 @@ class FESystem : public FiniteElement virtual FiniteElementDomination::Domination compare_for_face_domination (const FiniteElement &fe_other) const; + + /** + * For faces with non-standard + * face_orientation in 3D, the shape + * functions on faces have to be permutated + * in order to be combined with the correct + * dofs. This function returns a vector of + * integer values, that have to be added to + * the index of a shape function in order + * to result in the permutated index. Prior + * content of the vector @p shifts is + * erased. In 3D a vector of length @p + * dofs_per_quad is returned, in 2D and 1D + * there is no need for permutation and a + * vector of length 0 is returned. + */ + virtual + void + get_face_shape_function_shifts (std::vector &shifts) const; //@} /** diff --git a/deal.II/deal.II/include/fe/fe_values.h b/deal.II/deal.II/include/fe/fe_values.h index addb01b923..560e0d2226 100644 --- a/deal.II/deal.II/include/fe/fe_values.h +++ b/deal.II/deal.II/include/fe/fe_values.h @@ -263,7 +263,38 @@ class FEValuesData * non-zero components. */ std::vector shape_function_to_row_table; - + + /** + * Vector containing the permutation of + * shape functions necessary if the faces + * of a cell have the wrong + * face_orientation. This is computed only + * once. Actually, this does not contain + * the permutation itself but rather the + * shift of indices needed to calculate the + * permutation. + */ + std::vector shift_in_face_shape_functions; + + /** + * Vector containing the permutation of + * shape functions due to faces with + * non-standard face_orientation on a given + * cell (recomputed on each cell.) If all + * faces are oriented according to the + * standard, this is the identity mapping. + */ + std::vector permutated_shape_functions; + + /** + * Bool flag indicating the need to update + * the @p permutated_shape_functions vector + * on each cell. This is only necessary in + * 3d and if the finite element has + * shape_functions on the face. + */ + bool update_shape_function_permutation; + /** * Original update flags handed * to the constructor of @@ -626,7 +657,19 @@ class FEValuesBase : protected FEValuesData, const unsigned int point_no, const unsigned int component) const; - + /** + * If shape functions belong to a face in + * 3D, they have to be permutated, if the + * face has non-standard face + * orientation. This functuion takes an + * index of a shape function (on a standard + * cell) and returns the corresponding + * shape function on the real cell. + */ + unsigned int + shift_shape_function_index (const unsigned int i) const; + + //@} /// @name FunctionAccess Access to values of global finite element functions //@{ @@ -1623,6 +1666,12 @@ class FEValuesBase : protected FEValuesData, */ UpdateFlags compute_update_flags (const UpdateFlags update_flags) const; + /** + * Reinit the permutation of (face) shape + * functions to match the present cell. + */ + void reinit(); + private: /** * Copy constructor. Since @@ -2322,15 +2371,17 @@ double FEValuesBase::shape_value (const unsigned int i, const unsigned int j) const { + const unsigned int I=shift_shape_function_index(i); + Assert (this->update_flags & update_values, ExcAccessToUninitializedField()); - Assert (fe->is_primitive (i), - ExcShapeFunctionNotPrimitive(i)); + Assert (fe->is_primitive (I), + ExcShapeFunctionNotPrimitive(I)); // if the entire FE is primitive, // then we can take a short-cut: if (fe->is_primitive()) - return this->shape_values(i,j); + return this->shape_values(I,j); else // otherwise, use the mapping // between shape function numbers @@ -2341,7 +2392,7 @@ FEValuesBase::shape_value (const unsigned int i, // question to which vector // component the call of this // function refers - return this->shape_values(this->shape_function_to_row_table[i], j); + return this->shape_values(this->shape_function_to_row_table[I], j); } @@ -2353,6 +2404,8 @@ FEValuesBase::shape_value_component (const unsigned int i, const unsigned int j, const unsigned int component) const { + const unsigned int I=shift_shape_function_index(i); + Assert (this->update_flags & update_values, ExcAccessToUninitializedField()); Assert (component < fe->n_components(), @@ -2367,10 +2420,10 @@ FEValuesBase::shape_value_component (const unsigned int i, // system_to_component_table only // works if the shape function is // primitive): - if (fe->is_primitive(i)) + if (fe->is_primitive(I)) { - if (component == fe->system_to_component_index(i).first) - return this->shape_values(this->shape_function_to_row_table[i],j); + if (component == fe->system_to_component_index(I).first) + return this->shape_values(this->shape_function_to_row_table[I],j); else return 0; } @@ -2385,7 +2438,7 @@ FEValuesBase::shape_value_component (const unsigned int i, // whether the shape function // is non-zero at all within // this component: - if (fe->get_nonzero_components(i)[component] == false) + if (fe->get_nonzero_components(I)[component] == false) return 0.; // count how many non-zero @@ -2397,10 +2450,10 @@ FEValuesBase::shape_value_component (const unsigned int i, // shape function in the arrays // we index presently: const unsigned int - row = (this->shape_function_to_row_table[i] + row = (this->shape_function_to_row_table[I] + - std::count (fe->get_nonzero_components(i).begin(), - fe->get_nonzero_components(i).begin()+component, + std::count (fe->get_nonzero_components(I).begin(), + fe->get_nonzero_components(I).begin()+component, true)); return this->shape_values(row, j); }; @@ -2414,19 +2467,21 @@ const Tensor<1,dim> & FEValuesBase::shape_grad (const unsigned int i, const unsigned int j) const { + const unsigned int I=shift_shape_function_index(i); + Assert (this->update_flags & update_gradients, ExcAccessToUninitializedField()); - Assert (fe->is_primitive (i), - ExcShapeFunctionNotPrimitive(i)); + Assert (fe->is_primitive (I), + ExcShapeFunctionNotPrimitive(I)); Assert (ishape_gradients.size(), - ExcIndexRange (i, 0, this->shape_gradients.size())); + ExcIndexRange (I, 0, this->shape_gradients.size())); Assert (jshape_gradients[0].size(), ExcIndexRange (j, 0, this->shape_gradients[0].size())); // if the entire FE is primitive, // then we can take a short-cut: if (fe->is_primitive()) - return this->shape_gradients[i][j]; + return this->shape_gradients[I][j]; else // otherwise, use the mapping // between shape function numbers @@ -2437,7 +2492,7 @@ FEValuesBase::shape_grad (const unsigned int i, // question to which vector // component the call of this // function refers - return this->shape_gradients[this->shape_function_to_row_table[i]][j]; + return this->shape_gradients[this->shape_function_to_row_table[I]][j]; } @@ -2449,6 +2504,8 @@ FEValuesBase::shape_grad_component (const unsigned int i, const unsigned int j, const unsigned int component) const { + const unsigned int I=shift_shape_function_index(i); + Assert (this->update_flags & update_gradients, ExcAccessToUninitializedField()); Assert (component < fe->n_components(), @@ -2463,10 +2520,10 @@ FEValuesBase::shape_grad_component (const unsigned int i, // system_to_component_table only // works if the shape function is // primitive): - if (fe->is_primitive(i)) + if (fe->is_primitive(I)) { - if (component == fe->system_to_component_index(i).first) - return this->shape_gradients[this->shape_function_to_row_table[i]][j]; + if (component == fe->system_to_component_index(I).first) + return this->shape_gradients[this->shape_function_to_row_table[I]][j]; else return Tensor<1,dim>(); } @@ -2481,7 +2538,7 @@ FEValuesBase::shape_grad_component (const unsigned int i, // whether the shape function // is non-zero at all within // this component: - if (fe->get_nonzero_components(i)[component] == false) + if (fe->get_nonzero_components(I)[component] == false) return Tensor<1,dim>(); // count how many non-zero @@ -2493,10 +2550,10 @@ FEValuesBase::shape_grad_component (const unsigned int i, // shape function in the arrays // we index presently: const unsigned int - row = (this->shape_function_to_row_table[i] + row = (this->shape_function_to_row_table[I] + - std::count (fe->get_nonzero_components(i).begin(), - fe->get_nonzero_components(i).begin()+component, + std::count (fe->get_nonzero_components(I).begin(), + fe->get_nonzero_components(I).begin()+component, true)); return this->shape_gradients[row][j]; }; @@ -2510,19 +2567,21 @@ const Tensor<2,dim> & FEValuesBase::shape_2nd_derivative (const unsigned int i, const unsigned int j) const { + const unsigned int I=shift_shape_function_index(i); + Assert (this->update_flags & update_second_derivatives, ExcAccessToUninitializedField()); - Assert (fe->is_primitive (i), + Assert (fe->is_primitive (I), ExcShapeFunctionNotPrimitive(i)); - Assert (ishape_2nd_derivatives.size(), - ExcIndexRange (i, 0, this->shape_2nd_derivatives.size())); + Assert (Ishape_2nd_derivatives.size(), + ExcIndexRange (I, 0, this->shape_2nd_derivatives.size())); Assert (jshape_2nd_derivatives[0].size(), ExcIndexRange (j, 0, this->shape_2nd_derivatives[0].size())); // if the entire FE is primitive, // then we can take a short-cut: if (fe->is_primitive()) - return this->shape_2nd_derivatives[i][j]; + return this->shape_2nd_derivatives[I][j]; else // otherwise, use the mapping // between shape function numbers @@ -2533,7 +2592,7 @@ FEValuesBase::shape_2nd_derivative (const unsigned int i, // question to which vector // component the call of this // function refers - return this->shape_2nd_derivatives[this->shape_function_to_row_table[i]][j]; + return this->shape_2nd_derivatives[this->shape_function_to_row_table[I]][j]; } @@ -2545,6 +2604,8 @@ FEValuesBase::shape_2nd_derivative_component (const unsigned int i, const unsigned int j, const unsigned int component) const { + const unsigned int I=shift_shape_function_index(i); + Assert (this->update_flags & update_second_derivatives, ExcAccessToUninitializedField()); Assert (component < fe->n_components(), @@ -2559,10 +2620,10 @@ FEValuesBase::shape_2nd_derivative_component (const unsigned int i, // system_to_component_table only // works if the shape function is // primitive): - if (fe->is_primitive(i)) + if (fe->is_primitive(I)) { - if (component == fe->system_to_component_index(i).first) - return this->shape_2nd_derivatives[this->shape_function_to_row_table[i]][j]; + if (component == fe->system_to_component_index(I).first) + return this->shape_2nd_derivatives[this->shape_function_to_row_table[I]][j]; else return Tensor<2,dim>(); } @@ -2589,10 +2650,10 @@ FEValuesBase::shape_2nd_derivative_component (const unsigned int i, // shape function in the arrays // we index presently: const unsigned int - row = (this->shape_function_to_row_table[i] + row = (this->shape_function_to_row_table[I] + - std::count (fe->get_nonzero_components(i).begin(), - fe->get_nonzero_components(i).begin()+component, + std::count (fe->get_nonzero_components(I).begin(), + fe->get_nonzero_components(I).begin()+component, true)); return this->shape_2nd_derivatives[row][j]; }; @@ -2676,6 +2737,28 @@ FEValuesBase::JxW (const unsigned int i) const } + +template +inline +unsigned int +FEValuesBase::shift_shape_function_index (const unsigned int i) const +{ + // standard implementation for 1D and 2D + Assert(idofs_per_cell, ExcInternalError()); + return i; +} + +template <> +inline +unsigned int +FEValuesBase<3>::shift_shape_function_index (const unsigned int i) const +{ + Assert(idofs_per_cell, ExcInternalError()); + return this->permutated_shape_functions[i]; +} + + + /*------------------------ Inline functions: FEValues ----------------------------*/ @@ -2761,6 +2844,8 @@ FEFaceValuesBase::boundary_form (const unsigned int i) const return this->boundary_forms[i]; } + + #endif // DOXYGEN DEAL_II_NAMESPACE_CLOSE diff --git a/deal.II/deal.II/source/fe/fe.cc b/deal.II/deal.II/source/fe/fe.cc index 199a8b9ab3..d08e3ff599 100644 --- a/deal.II/deal.II/source/fe/fe.cc +++ b/deal.II/deal.II/source/fe/fe.cc @@ -333,6 +333,30 @@ FiniteElement::component_to_block_index (const unsigned int index) const } +#if deal_II_dimension < 3 + +template +void +FiniteElement::get_face_shape_function_shifts (std::vector &shifts) const +{ + // general template for 1D and 2D, return an + // empty vector + shifts.clear(); +} + +#else + +template <> +void +FiniteElement<3>::get_face_shape_function_shifts (std::vector &shifts) const +{ + shifts.clear(); + shifts.resize(this->dofs_per_quad,0); +} + +#endif + + template bool diff --git a/deal.II/deal.II/source/fe/fe_q.cc b/deal.II/deal.II/source/fe/fe_q.cc index ceb3474d86..b584801c38 100644 --- a/deal.II/deal.II/source/fe/fe_q.cc +++ b/deal.II/deal.II/source/fe/fe_q.cc @@ -710,6 +710,37 @@ compare_for_face_domination (const FiniteElement &fe_other) const +template +void +FE_Q::get_face_shape_function_shifts (std::vector &shifts) const +{ + // general template for 1D and 2D, return an + // empty vector + shifts.clear(); +} + + + +#if deal_II_dimension == 3 + +template <> +void +FE_Q<3>::get_face_shape_function_shifts (std::vector &shifts) const +{ + shifts.resize(this->dofs_per_quad); + + unsigned int points=this->degree-1; + Assert(points*points==this->dofs_per_quad, ExcInternalError()); + + for (unsigned int local=0; localdofs_per_quad; ++local) + // face support points are in lexicographic + // ordering with x running fastest. invert + // that (y running fastest) + shifts[local] = (local%points)*points + local/points - local; +} + +#endif + //--------------------------------------------------------------------------- // Auxiliary functions //--------------------------------------------------------------------------- diff --git a/deal.II/deal.II/source/fe/fe_system.cc b/deal.II/deal.II/source/fe/fe_system.cc index 74db54ecd3..48f2e4b7a8 100644 --- a/deal.II/deal.II/source/fe/fe_system.cc +++ b/deal.II/deal.II/source/fe/fe_system.cc @@ -2926,6 +2926,40 @@ FESystem::unit_face_support_point (const unsigned index) const +template +void +FESystem::get_face_shape_function_shifts (std::vector &shifts) const +{ + // general template for 1D and 2D, return an + // empty vector + shifts.clear(); +} + + + +#if deal_II_dimension == 3 + +template <> +void +FESystem<3>::get_face_shape_function_shifts (std::vector &shifts) const +{ + shifts.clear(); + std::vector temp; + // to obtain the shifts for this composed + // element, concatenate the shift vectors of + // the base elements + for (unsigned int b=0; bbase_element(b).get_face_shape_function_shifts(temp); + for (unsigned int c=0; cdofs_per_quad, ExcInternalError()); +} + +#endif + + template unsigned int diff --git a/deal.II/deal.II/source/fe/fe_values.cc b/deal.II/deal.II/source/fe/fe_values.cc index 84aeac4bf5..879738a547 100644 --- a/deal.II/deal.II/source/fe/fe_values.cc +++ b/deal.II/deal.II/source/fe/fe_values.cc @@ -316,6 +316,34 @@ FEValuesData::initialize (const unsigned int n_quadrature_points, if (flags & update_cell_JxW_values) this->cell_JxW_values.resize(n_quadrature_points); + + // initialize the permutation fields, if they + // are needed + if (dim==3) + { + permutated_shape_functions.resize(fe.dofs_per_cell); + // initialize cell permutation mapping + // with identity + for (unsigned int i=0; ipermutated_shape_functions[i]=i; + + if (fe.dofs_per_quad>0 && + (flags & update_values || + flags & update_gradients || + flags & update_second_derivatives)) + { + // ask the fe to fill the vector of + // shifts + fe.get_face_shape_function_shifts(shift_in_face_shape_functions); + Assert (shift_in_face_shape_functions.size()==fe.dofs_per_quad, + ExcInternalError()); + update_shape_function_permutation=true; + } + else + update_shape_function_permutation=false; + } + else + update_shape_function_permutation=false; } @@ -365,6 +393,43 @@ FEValuesBase::~FEValuesBase () } +#if deal_II_dimension <3 +template +void FEValuesBase::reinit () +{ + // do nothing in 1D and 2D +} +#else + +template <> +void FEValuesBase<3>::reinit () +{ + // in 3D reinit the permutation of face shape + // functions, if necessary + if (update_shape_function_permutation) + { + Assert (this->shift_in_face_shape_functions.size()==this->fe->dofs_per_quad, + ExcInternalError()); + Triangulation<3>::cell_iterator thiscell = *(this->present_cell); + unsigned int offset=fe->first_quad_index; + for (unsigned int face_no=0; face_no::faces_per_cell; ++face_no) + { + if (thiscell->face_orientation(face_no)) + for (unsigned int i=0; idofs_per_quad; ++i) + this->permutated_shape_functions[offset+i]=offset+i; + else + for (unsigned int i=0; idofs_per_quad; ++i) + this->permutated_shape_functions[offset+i]=offset+i+this->shift_in_face_shape_functions[i]; + + offset+=this->fe->dofs_per_quad; + } + Assert (offset-fe->first_quad_index==GeometryInfo<3>::faces_per_cell*this->fe->dofs_per_quad, + ExcInternalError()); + } +} +#endif + + template template @@ -1011,12 +1076,15 @@ FEValuesBase::memory_consumption () const return (MemoryConsumption::memory_consumption (this->shape_values) + MemoryConsumption::memory_consumption (this->shape_gradients) + MemoryConsumption::memory_consumption (this->shape_2nd_derivatives) + + MemoryConsumption::memory_consumption (this->permutated_shape_functions) + + MemoryConsumption::memory_consumption (this->shift_in_face_shape_functions) + MemoryConsumption::memory_consumption (this->JxW_values) + MemoryConsumption::memory_consumption (this->quadrature_points) + MemoryConsumption::memory_consumption (this->normal_vectors) + MemoryConsumption::memory_consumption (this->boundary_forms) + MemoryConsumption::memory_consumption (this->cell_JxW_values) + sizeof(this->update_flags) + + sizeof(this->update_shape_function_permutation) + MemoryConsumption::memory_consumption (n_quadrature_points) + MemoryConsumption::memory_consumption (dofs_per_cell) + MemoryConsumption::memory_consumption (mapping) + @@ -1260,6 +1328,8 @@ void FEValues::do_reinit () this->fe_data->clear_first_cell (); this->mapping_data->clear_first_cell (); + + FEValuesBase::reinit(); } @@ -1537,6 +1607,7 @@ void FEFaceValues::do_reinit (const unsigned int face_no) this->fe_data->clear_first_cell (); this->mapping_data->clear_first_cell (); + FEValuesBase::reinit(); } @@ -1789,6 +1860,7 @@ void FESubfaceValues::do_reinit (const unsigned int face_no, this->fe_data->clear_first_cell (); this->mapping_data->clear_first_cell (); + FEValuesBase::reinit(); } diff --git a/deal.II/doc/news/changes.html b/deal.II/doc/news/changes.html index eee625f1a2..77293dbeaa 100644 --- a/deal.II/doc/news/changes.html +++ b/deal.II/doc/news/changes.html @@ -959,6 +959,16 @@ inconvenience this causes.

deal.II

    +
  1. Fixed: On faces with wrong face_orientation the shape + functions have to reordered in order to be combined with the correct + dofs. This is only relevant for continuous elements in 3D. At least for + FE_Q and systems of FE_Q this works now, for other finite elements the + reordering vector still has to be implemented. +
    + (Tobias Leicht, 2007/01/17) +

    +
  2. Fixed: The Triangulation::execute_coarsening_and_refinement function has to -- 2.39.5