From: grahambenharper Date: Sun, 7 Jul 2019 03:30:22 +0000 (-0600) Subject: Add functionality to support multiple mappings per element in FE_PolyTensor X-Git-Tag: v9.2.0-rc1~1381^2~4 X-Git-Url: https://gitweb.dealii.org/cgi-bin/gitweb.cgi?a=commitdiff_plain;h=00fa59797264ad174588ea35ed439252b5ea0197;p=dealii.git Add functionality to support multiple mappings per element in FE_PolyTensor --- diff --git a/include/deal.II/fe/fe_dg_vector.templates.h b/include/deal.II/fe/fe_dg_vector.templates.h index 8aa11a5a29..a0893d2e47 100644 --- a/include/deal.II/fe/fe_dg_vector.templates.h +++ b/include/deal.II/fe/fe_dg_vector.templates.h @@ -43,7 +43,7 @@ FE_DGVector::FE_DGVector(const unsigned int deg, std::vector(PolynomialType::compute_n_pols(deg), ComponentMask(dim, true))) { - this->mapping_type = map; + this->mapping_type = {map}; const unsigned int polynomial_degree = this->tensor_degree(); QGauss quadrature(polynomial_degree + 1); diff --git a/include/deal.II/fe/fe_poly_tensor.h b/include/deal.II/fe/fe_poly_tensor.h index ecff871dc7..c755a6f2a6 100644 --- a/include/deal.II/fe/fe_poly_tensor.h +++ b/include/deal.II/fe/fe_poly_tensor.h @@ -207,9 +207,12 @@ public: protected: /** * The mapping type to be used to map shape functions from the reference - * cell to the mesh cell. + * cell to the mesh cell. If this vector is length one, the same mapping + * will be applied to all shape functions. If the vector size is equal to + * the finite element dofs per cell, then each shape function will be mapped + * according to the corresponding entry in the vector. */ - MappingType mapping_type; + std::vector mapping_type; /* NOTE: The following function has its definition inlined into the class @@ -247,11 +250,30 @@ protected: // initialize fields only if really // necessary. otherwise, don't // allocate memory + + const bool update_transformed_shape_values = + (std::find_if(this->mapping_type.begin(), + this->mapping_type.end(), + [](const MappingType t) { return t != mapping_none; }) != + this->mapping_type.end()); + + const bool update_transformed_shape_grads = + (std::find_if(this->mapping_type.begin(), + this->mapping_type.end(), + [](const MappingType t) { + return (t == mapping_raviart_thomas || + t == mapping_piola || t == mapping_nedelec || + t == mapping_contravariant); + }) != this->mapping_type.end()); + + const bool update_transformed_shape_hessian_tensors = + update_transformed_shape_values; + if (update_flags & update_values) { values.resize(this->dofs_per_cell); data.shape_values.reinit(this->dofs_per_cell, n_q_points); - if (mapping_type != mapping_none) + if (update_transformed_shape_values) data.transformed_shape_values.resize(n_q_points); } @@ -261,10 +283,7 @@ protected: data.shape_grads.reinit(this->dofs_per_cell, n_q_points); data.transformed_shape_grads.resize(n_q_points); - if ((mapping_type == mapping_raviart_thomas) || - (mapping_type == mapping_piola) || - (mapping_type == mapping_nedelec) || - (mapping_type == mapping_contravariant)) + if (update_transformed_shape_grads) data.untransformed_shape_grads.resize(n_q_points); } @@ -273,7 +292,7 @@ protected: grad_grads.resize(this->dofs_per_cell); data.shape_grad_grads.reinit(this->dofs_per_cell, n_q_points); data.transformed_shape_hessians.resize(n_q_points); - if (mapping_type != mapping_none) + if (update_transformed_shape_hessian_tensors) data.untransformed_shape_hessian_tensors.resize(n_q_points); } diff --git a/source/fe/fe_abf.cc b/source/fe/fe_abf.cc index 9efdc6e053..229344046f 100644 --- a/source/fe/fe_abf.cc +++ b/source/fe/fe_abf.cc @@ -60,7 +60,7 @@ FE_ABF::FE_ABF(const unsigned int deg) Assert(dim >= 2, ExcImpossibleInDim(dim)); const unsigned int n_dofs = this->dofs_per_cell; - this->mapping_type = mapping_raviart_thomas; + this->mapping_type = {mapping_raviart_thomas}; // First, initialize the // generalized support points and // quadrature weights, since they diff --git a/source/fe/fe_bdm.cc b/source/fe/fe_bdm.cc index 2b5335d940..0dfd5374f6 100644 --- a/source/fe/fe_bdm.cc +++ b/source/fe/fe_bdm.cc @@ -56,7 +56,7 @@ FE_BDM::FE_BDM(const unsigned int deg) const unsigned int n_dofs = this->dofs_per_cell; - this->mapping_type = mapping_bdm; + this->mapping_type = {mapping_bdm}; // These must be done first, since // they change the evaluation of // basis functions diff --git a/source/fe/fe_dg_vector.cc b/source/fe/fe_dg_vector.cc index 229c009ddb..9409a368cb 100644 --- a/source/fe/fe_dg_vector.cc +++ b/source/fe/fe_dg_vector.cc @@ -27,7 +27,7 @@ DEAL_II_NAMESPACE_OPEN template FE_DGNedelec::FE_DGNedelec(const unsigned int p) - : FE_DGVector, dim, spacedim>(p, mapping_nedelec) + : FE_DGVector, dim, spacedim>(p, {mapping_nedelec}) {} @@ -52,7 +52,7 @@ template FE_DGRaviartThomas::FE_DGRaviartThomas(const unsigned int p) : FE_DGVector, dim, spacedim>( p, - mapping_raviart_thomas) + {mapping_raviart_thomas}) {} @@ -77,7 +77,7 @@ FE_DGRaviartThomas::get_name() const template FE_DGBDM::FE_DGBDM(const unsigned int p) - : FE_DGVector, dim, spacedim>(p, mapping_bdm) + : FE_DGVector, dim, spacedim>(p, {mapping_bdm}) {} diff --git a/source/fe/fe_nedelec.cc b/source/fe/fe_nedelec.cc index d12afe8a95..ad175ac708 100644 --- a/source/fe/fe_nedelec.cc +++ b/source/fe/fe_nedelec.cc @@ -87,7 +87,7 @@ FE_Nedelec::FE_Nedelec(const unsigned int order) const unsigned int n_dofs = this->dofs_per_cell; - this->mapping_type = mapping_nedelec; + this->mapping_type = {mapping_nedelec}; // First, initialize the // generalized support points and // quadrature weights, since they diff --git a/source/fe/fe_poly_tensor.cc b/source/fe/fe_poly_tensor.cc index 6b72947c55..bc1d0844a3 100644 --- a/source/fe/fe_poly_tensor.cc +++ b/source/fe/fe_poly_tensor.cc @@ -49,30 +49,29 @@ namespace internal * implements an algorithm that determines those DoFs that need this * sign change for a given cell. */ + template void get_face_sign_change_rt(const dealii::Triangulation<1>::cell_iterator &, - const unsigned int, - std::vector &face_sign) + const FiniteElement<1, spacedim> &, + const std::vector &, + std::vector &) { // nothing to do in 1d - std::fill(face_sign.begin(), face_sign.end(), 1.0); } + // template void get_face_sign_change_rt( const dealii::Triangulation<2>::cell_iterator &cell, - const unsigned int dofs_per_face, + const FiniteElement<2, 2> & fe, + const std::vector & mapping_type, std::vector & face_sign) { const unsigned int dim = 2; const unsigned int spacedim = 2; - // Default is no sign - // change. I.e. multiply by one. - std::fill(face_sign.begin(), face_sign.end(), 1.0); - for (unsigned int f = GeometryInfo::faces_per_cell / 2; f < GeometryInfo::faces_per_cell; ++f) @@ -84,14 +83,24 @@ namespace internal const unsigned int nn = cell->neighbor_face_no(f); if (nn < GeometryInfo::faces_per_cell / 2) - for (unsigned int j = 0; j < dofs_per_face; ++j) + for (unsigned int j = 0; j < fe.dofs_per_face; ++j) { - Assert(f * dofs_per_face + j < face_sign.size(), + const unsigned int cell_j = fe.face_to_cell_index(j, f); + // something in FiniteElement or FiniteElementData base + // classes + + Assert(f * fe.dofs_per_face + j < face_sign.size(), + ExcInternalError()); + Assert(mapping_type.size() == 1 || + cell_j < mapping_type.size(), ExcInternalError()); // TODO: This is probably only going to work for those // elements for which all dofs are face dofs - face_sign[f * dofs_per_face + j] = -1.0; + if ((mapping_type.size() > 1 ? + mapping_type[cell_j] : + mapping_type[0]) == mapping_raviart_thomas) + face_sign[f * fe.dofs_per_face + j] = -1.0; } } } @@ -99,51 +108,57 @@ namespace internal + template void get_face_sign_change_rt( const dealii::Triangulation<3>::cell_iterator & /*cell*/, - const unsigned int /*dofs_per_face*/, - std::vector &face_sign) + const FiniteElement<3, spacedim> & /*fe*/, + const std::vector & /*mapping_type*/, + std::vector & /*face_sign*/) { - std::fill(face_sign.begin(), face_sign.end(), 1.0); // TODO: think about what it would take here } + + + template void get_face_sign_change_nedelec( const dealii::Triangulation<1>::cell_iterator & /*cell*/, - const unsigned int /*dofs_per_face*/, - std::vector &face_sign) + const FiniteElement<1, spacedim> & /*fe*/, + const std::vector & /*mapping_type*/, + std::vector & /*face_sign*/) { // nothing to do in 1d - std::fill(face_sign.begin(), face_sign.end(), 1.0); } + template void get_face_sign_change_nedelec( const dealii::Triangulation<2>::cell_iterator & /*cell*/, - const unsigned int /*dofs_per_face*/, - std::vector &face_sign) + const FiniteElement<2, spacedim> & /*fe*/, + const std::vector & /*mapping_type*/, + std::vector & /*face_sign*/) { - std::fill(face_sign.begin(), face_sign.end(), 1.0); // TODO: think about what it would take here } + template void get_face_sign_change_nedelec( const dealii::Triangulation<3>::cell_iterator &cell, - const unsigned int /*dofs_per_face*/, - std::vector &face_sign) + const FiniteElement<3, spacedim> & /*fe*/, + const std::vector &mapping_type, + std::vector & face_sign) { const unsigned int dim = 3; - std::fill(face_sign.begin(), face_sign.end(), 1.0); // TODO: This is probably only going to work for those elements for // which all dofs are face dofs for (unsigned int l = 0; l < GeometryInfo::lines_per_cell; ++l) - if (!(cell->line_orientation(l))) + if (!(cell->line_orientation(l)) & mapping_type[0] == mapping_nedelec) face_sign[l] = -1.0; } } // namespace @@ -161,7 +176,7 @@ FE_PolyTensor::FE_PolyTensor( : FiniteElement(fe_data, restriction_is_additive_flags, nonzero_components) - , mapping_type(MappingType::mapping_none) + , mapping_type({MappingType::mapping_none}) , poly_space(PolynomialType(degree)) { cached_point(0) = -1; @@ -364,18 +379,23 @@ FE_PolyTensor::fill_fe_values( // between two faces. std::fill(fe_data.sign_change.begin(), fe_data.sign_change.end(), 1.0); - if (mapping_type == mapping_raviart_thomas) - internal::FE_PolyTensor::get_face_sign_change_rt(cell, - this->dofs_per_face, - fe_data.sign_change); - else if (mapping_type == mapping_nedelec) - internal::FE_PolyTensor::get_face_sign_change_nedelec(cell, - this->dofs_per_face, - fe_data.sign_change); + internal::FE_PolyTensor::get_face_sign_change_rt(cell, + *this, + this->mapping_type, + fe_data.sign_change); + + internal::FE_PolyTensor::get_face_sign_change_nedelec(cell, + *this, + this->mapping_type, + fe_data.sign_change); for (unsigned int i = 0; i < this->dofs_per_cell; ++i) { + MappingType mapping_type = + (this->mapping_type.size() > 1 ? this->mapping_type[i] : + this->mapping_type[0]); + const unsigned int first = output_data.shape_function_to_row_table[i * this->n_components() + this->get_nonzero_components(i) @@ -940,18 +960,22 @@ FE_PolyTensor::fill_fe_face_values( // on the neighborhood between two faces. std::fill(fe_data.sign_change.begin(), fe_data.sign_change.end(), 1.0); - if (mapping_type == mapping_raviart_thomas) - internal::FE_PolyTensor::get_face_sign_change_rt(cell, - this->dofs_per_face, - fe_data.sign_change); + internal::FE_PolyTensor::get_face_sign_change_rt(cell, + *this, + this->mapping_type, + fe_data.sign_change); - else if (mapping_type == mapping_nedelec) - internal::FE_PolyTensor::get_face_sign_change_nedelec(cell, - this->dofs_per_face, - fe_data.sign_change); + internal::FE_PolyTensor::get_face_sign_change_nedelec(cell, + *this, + this->mapping_type, + fe_data.sign_change); for (unsigned int i = 0; i < this->dofs_per_cell; ++i) { + const MappingType mapping_type = + (this->mapping_type.size() > 1 ? this->mapping_type[i] : + this->mapping_type[0]); + const unsigned int first = output_data.shape_function_to_row_table[i * this->n_components() + this->get_nonzero_components(i) @@ -1568,18 +1592,22 @@ FE_PolyTensor::fill_fe_subface_values( // on the neighborhood between two faces. std::fill(fe_data.sign_change.begin(), fe_data.sign_change.end(), 1.0); - if (mapping_type == mapping_raviart_thomas) - internal::FE_PolyTensor::get_face_sign_change_rt(cell, - this->dofs_per_face, - fe_data.sign_change); + internal::FE_PolyTensor::get_face_sign_change_rt(cell, + *this, + this->mapping_type, + fe_data.sign_change); - else if (mapping_type == mapping_nedelec) - internal::FE_PolyTensor::get_face_sign_change_nedelec(cell, - this->dofs_per_face, - fe_data.sign_change); + internal::FE_PolyTensor::get_face_sign_change_nedelec(cell, + *this, + this->mapping_type, + fe_data.sign_change); for (unsigned int i = 0; i < this->dofs_per_cell; ++i) { + MappingType mapping_type = + (this->mapping_type.size() > 1 ? this->mapping_type[i] : + this->mapping_type[0]); + const unsigned int first = output_data.shape_function_to_row_table[i * this->n_components() + this->get_nonzero_components(i) @@ -2140,88 +2168,95 @@ FE_PolyTensor::requires_update_flags( { UpdateFlags out = update_default; - switch (mapping_type) + for (unsigned int i = 0; i < this->dofs_per_cell; ++i) { - case mapping_none: - { - if (flags & update_values) - out |= update_values; - - if (flags & update_gradients) - out |= update_gradients | update_values | - update_jacobian_pushed_forward_grads; - - if (flags & update_hessians) - out |= update_hessians | update_values | update_gradients | - update_jacobian_pushed_forward_grads | - update_jacobian_pushed_forward_2nd_derivatives; - break; - } + MappingType mapping_type = + (this->mapping_type.size() > 1 ? this->mapping_type[i] : + this->mapping_type[0]); - case mapping_raviart_thomas: - case mapping_piola: + switch (mapping_type) { - if (flags & update_values) - out |= update_values | update_piola; - - if (flags & update_gradients) - out |= update_gradients | update_values | update_piola | - update_jacobian_pushed_forward_grads | - update_covariant_transformation | - update_contravariant_transformation; - - if (flags & update_hessians) - out |= update_hessians | update_piola | update_values | - update_gradients | update_jacobian_pushed_forward_grads | - update_jacobian_pushed_forward_2nd_derivatives | - update_covariant_transformation; - - break; - } + case mapping_none: + { + if (flags & update_values) + out |= update_values; + + if (flags & update_gradients) + out |= update_gradients | update_values | + update_jacobian_pushed_forward_grads; + + if (flags & update_hessians) + out |= update_hessians | update_values | update_gradients | + update_jacobian_pushed_forward_grads | + update_jacobian_pushed_forward_2nd_derivatives; + break; + } + case mapping_raviart_thomas: + case mapping_piola: + { + if (flags & update_values) + out |= update_values | update_piola; + + if (flags & update_gradients) + out |= update_gradients | update_values | update_piola | + update_jacobian_pushed_forward_grads | + update_covariant_transformation | + update_contravariant_transformation; + + if (flags & update_hessians) + out |= update_hessians | update_piola | update_values | + update_gradients | update_jacobian_pushed_forward_grads | + update_jacobian_pushed_forward_2nd_derivatives | + update_covariant_transformation; + + break; + } - case mapping_contravariant: - { - if (flags & update_values) - out |= update_values | update_piola; - - if (flags & update_gradients) - out |= update_gradients | update_values | - update_jacobian_pushed_forward_grads | - update_covariant_transformation | - update_contravariant_transformation; - - if (flags & update_hessians) - out |= update_hessians | update_piola | update_values | - update_gradients | update_jacobian_pushed_forward_grads | - update_jacobian_pushed_forward_2nd_derivatives | - update_covariant_transformation; - - break; - } - case mapping_nedelec: - case mapping_covariant: - { - if (flags & update_values) - out |= update_values | update_covariant_transformation; + case mapping_contravariant: + { + if (flags & update_values) + out |= update_values | update_piola; + + if (flags & update_gradients) + out |= update_gradients | update_values | + update_jacobian_pushed_forward_grads | + update_covariant_transformation | + update_contravariant_transformation; + + if (flags & update_hessians) + out |= update_hessians | update_piola | update_values | + update_gradients | update_jacobian_pushed_forward_grads | + update_jacobian_pushed_forward_2nd_derivatives | + update_covariant_transformation; + + break; + } - if (flags & update_gradients) - out |= update_gradients | update_values | - update_jacobian_pushed_forward_grads | - update_covariant_transformation; + case mapping_nedelec: + case mapping_covariant: + { + if (flags & update_values) + out |= update_values | update_covariant_transformation; - if (flags & update_hessians) - out |= update_hessians | update_values | update_gradients | - update_jacobian_pushed_forward_grads | - update_jacobian_pushed_forward_2nd_derivatives | - update_covariant_transformation; + if (flags & update_gradients) + out |= update_gradients | update_values | + update_jacobian_pushed_forward_grads | + update_covariant_transformation; - break; - } + if (flags & update_hessians) + out |= update_hessians | update_values | update_gradients | + update_jacobian_pushed_forward_grads | + update_jacobian_pushed_forward_2nd_derivatives | + update_covariant_transformation; - default: - { - Assert(false, ExcNotImplemented()); + break; + } + + default: + { + Assert(false, ExcNotImplemented()); + } } } diff --git a/source/fe/fe_raviart_thomas.cc b/source/fe/fe_raviart_thomas.cc index 779623ee62..0aead6d54e 100644 --- a/source/fe/fe_raviart_thomas.cc +++ b/source/fe/fe_raviart_thomas.cc @@ -60,7 +60,7 @@ FE_RaviartThomas::FE_RaviartThomas(const unsigned int deg) Assert(dim >= 2, ExcImpossibleInDim(dim)); const unsigned int n_dofs = this->dofs_per_cell; - this->mapping_type = mapping_raviart_thomas; + this->mapping_type = {mapping_raviart_thomas}; // First, initialize the // generalized support points and // quadrature weights, since they diff --git a/source/fe/fe_raviart_thomas_nodal.cc b/source/fe/fe_raviart_thomas_nodal.cc index 7368158080..a7f40dfc98 100644 --- a/source/fe/fe_raviart_thomas_nodal.cc +++ b/source/fe/fe_raviart_thomas_nodal.cc @@ -52,7 +52,7 @@ FE_RaviartThomasNodal::FE_RaviartThomasNodal(const unsigned int deg) Assert(dim >= 2, ExcImpossibleInDim(dim)); const unsigned int n_dofs = this->dofs_per_cell; - this->mapping_type = mapping_raviart_thomas; + this->mapping_type = {mapping_raviart_thomas}; // First, initialize the // generalized support points and // quadrature weights, since they diff --git a/source/fe/fe_rt_bubbles.cc b/source/fe/fe_rt_bubbles.cc index 240e978698..a5ab697f27 100644 --- a/source/fe/fe_rt_bubbles.cc +++ b/source/fe/fe_rt_bubbles.cc @@ -55,7 +55,7 @@ FE_RT_Bubbles::FE_RT_Bubbles(const unsigned int deg) "Lowest order RT_Bubbles element is degree 1, but you requested for degree 0")); const unsigned int n_dofs = this->dofs_per_cell; - this->mapping_type = mapping_raviart_thomas; + this->mapping_type = {mapping_raviart_thomas}; // Initialize support points and quadrature weights initialize_support_points(deg); // Compute the inverse node matrix to get