From 7ce0f37f354e389f9813de7e322b01a8c2df86ba Mon Sep 17 00:00:00 2001 From: Nicola Giuliani Date: Thu, 24 Mar 2016 23:32:10 +0100 Subject: [PATCH] seems to work --- source/fe/fe_tools.cc | 1050 +++++++++++++++++++++-------------------- 1 file changed, 529 insertions(+), 521 deletions(-) diff --git a/source/fe/fe_tools.cc b/source/fe/fe_tools.cc index b1af515850..9fc4c1a710 100644 --- a/source/fe/fe_tools.cc +++ b/source/fe/fe_tools.cc @@ -97,11 +97,24 @@ namespace FETools return new FE_DGQArbitraryNodes<1>(quad); } template <> + FiniteElement<1, 2> * + FEFactory >::get (const Quadrature<1> &quad) const + { + return new FE_DGQArbitraryNodes<1, 2>(quad); + } + template <> FiniteElement<2, 2> * FEFactory >::get (const Quadrature<1> &quad) const { return new FE_DGQArbitraryNodes<2>(quad); } + template <> + FiniteElement<2, 3> * + FEFactory >::get (const Quadrature<1> &quad) const + { + return new FE_DGQArbitraryNodes<2, 3>(quad); + } + template <> FiniteElement<3, 3> * FEFactory >::get (const Quadrature<1> &quad) const @@ -1511,617 +1524,612 @@ namespace FETools const FEFactoryBase *fef=dynamic_cast*>(ptr); return fef->get(QGaussLobatto<1>(tmp.first)); } + else if (quadrature_name.compare("QGauss") == 0) + { + const std::pair tmp + = Utilities::get_integer_at_position (name, 0); + // delete "))" + name.erase(0, tmp.second+2); + const Subscriptor *ptr = fe_name_map.find(name_part)->second.get(); + const FEFactoryBase *fef=dynamic_cast*>(ptr); + return fef->get(QGauss<1>(tmp.first)); + } else { - unsigned int position = name.find('('); - const std::string quadrature_name(name, 0, position); - name.erase(0,position+1); - if (quadrature_name.compare("QGauss") == 0) - { - const std::pair tmp - = Utilities::get_integer_at_position (name, 0); - // delete "))" - name.erase(0, tmp.second+2); - const Subscriptor *ptr = fe_name_map.find(name_part)->second.get(); - const FEFactoryBase *fef=dynamic_cast*>(ptr); - return fef->get(QGauss<1>(tmp.first)); - } - else - { - AssertThrow (false,ExcNotImplemented()); - } + AssertThrow (false,ExcNotImplemented()); } } + } - // hm, if we have come thus far, we - // didn't know what to do with the - // string we got. so do as the docs - // say: raise an exception - AssertThrow (false, ExcInvalidFEName(name)); + // hm, if we have come thus far, we + // didn't know what to do with the + // string we got. so do as the docs + // say: raise an exception + AssertThrow (false, ExcInvalidFEName(name)); - // make some compilers happy that - // do not realize that we can't get - // here after throwing - return 0; - } + // make some compilers happy that + // do not realize that we can't get + // here after throwing + return 0; + } - template - FiniteElement *get_fe_from_name (std::string &name) - { - return get_fe_from_name_ext (name, fe_name_map[dim][spacedim]); - } + template + FiniteElement *get_fe_from_name (std::string &name) + { + return get_fe_from_name_ext (name, fe_name_map[dim][spacedim]); } } + } - template - FiniteElement * - get_fe_by_name (const std::string ¶meter_name) - { - std::string name = Utilities::trim(parameter_name); - std::size_t index = 1; - // remove spaces that are not between two word (things that match the - // regular expression [A-Za-z0-9_]) characters. - while (2 < name.size() && index < name.size() - 1) - { - if (name[index] == ' ' && - (!(std::isalnum(name[index - 1]) || name[index - 1] == '_') || - !(std::isalnum(name[index + 1]) || name[index + 1] == '_'))) - { - name.erase(index, 1); - } - else - { - ++index; - } - } + template + FiniteElement * + get_fe_by_name (const std::string ¶meter_name) + { + std::string name = Utilities::trim(parameter_name); + std::size_t index = 1; + // remove spaces that are not between two word (things that match the + // regular expression [A-Za-z0-9_]) characters. + while (2 < name.size() && index < name.size() - 1) + { + if (name[index] == ' ' && + (!(std::isalnum(name[index - 1]) || name[index - 1] == '_') || + !(std::isalnum(name[index + 1]) || name[index + 1] == '_'))) + { + name.erase(index, 1); + } + else + { + ++index; + } + } - // Create a version of the name - // string where all template - // parameters are eliminated. - for (unsigned int pos1 = name.find('<'); - pos1 < name.size(); - pos1 = name.find('<')) - { + // Create a version of the name + // string where all template + // parameters are eliminated. + for (unsigned int pos1 = name.find('<'); + pos1 < name.size(); + pos1 = name.find('<')) + { - const unsigned int pos2 = name.find('>'); - // If there is only a single - // character between those two, - // it should be 'd' or the number - // representing the dimension. - if (pos2-pos1 == 2) - { - const char dimchar = '0' + dim; - (void)dimchar; - if (name.at(pos1+1) != 'd') - Assert (name.at(pos1+1) == dimchar, - ExcInvalidFEDimension(name.at(pos1+1), dim)); - } - else - Assert(pos2-pos1 == 4, ExcInvalidFEName(name)); + const unsigned int pos2 = name.find('>'); + // If there is only a single + // character between those two, + // it should be 'd' or the number + // representing the dimension. + if (pos2-pos1 == 2) + { + const char dimchar = '0' + dim; + (void)dimchar; + if (name.at(pos1+1) != 'd') + Assert (name.at(pos1+1) == dimchar, + ExcInvalidFEDimension(name.at(pos1+1), dim)); + } + else + Assert(pos2-pos1 == 4, ExcInvalidFEName(name)); - // If pos1==pos2, then we are - // probably at the end of the - // string - if (pos2 != pos1) - name.erase(pos1, pos2-pos1+1); - } - // Replace all occurrences of "^dim" - // by "^d" to be handled by the - // next loop - for (unsigned int pos = name.find("^dim"); - pos < name.size(); - pos = name.find("^dim")) - name.erase(pos+2, 2); - - // Replace all occurrences of "^d" - // by using the actual dimension - for (unsigned int pos = name.find("^d"); - pos < name.size(); - pos = name.find("^d")) - name.at(pos+1) = '0' + dim; - - try - { - FiniteElement *fe = internal::get_fe_from_name (name); - - // Make sure the auxiliary function - // ate up all characters of the name. - AssertThrow (name.size() == 0, - ExcInvalidFEName(parameter_name - + std::string(" extra characters after " - "end of name"))); - return fe; - } - catch (const std::string &errline) - { - AssertThrow(false, ExcInvalidFEName(parameter_name - + std::string(" at ") - + errline)); - return 0; - } - } + // If pos1==pos2, then we are + // probably at the end of the + // string + if (pos2 != pos1) + name.erase(pos1, pos2-pos1+1); + } + // Replace all occurrences of "^dim" + // by "^d" to be handled by the + // next loop + for (unsigned int pos = name.find("^dim"); + pos < name.size(); + pos = name.find("^dim")) + name.erase(pos+2, 2); + + // Replace all occurrences of "^d" + // by using the actual dimension + for (unsigned int pos = name.find("^d"); + pos < name.size(); + pos = name.find("^d")) + name.at(pos+1) = '0' + dim; + try + { + FiniteElement *fe = internal::get_fe_from_name (name); + + // Make sure the auxiliary function + // ate up all characters of the name. + AssertThrow (name.size() == 0, + ExcInvalidFEName(parameter_name + + std::string(" extra characters after " + "end of name"))); + return fe; + } + catch (const std::string &errline) + { + AssertThrow(false, ExcInvalidFEName(parameter_name + + std::string(" at ") + + errline)); + return 0; + } + } - template - FiniteElement * - get_fe_from_name (const std::string ¶meter_name) - { - return get_fe_by_name (parameter_name); - } + template + FiniteElement * + get_fe_from_name (const std::string ¶meter_name) + { + return get_fe_by_name (parameter_name); + } - template - void - compute_projection_from_quadrature_points_matrix (const FiniteElement &fe, - const Quadrature &lhs_quadrature, - const Quadrature &rhs_quadrature, - FullMatrix &X) - { - Assert (fe.n_components() == 1, ExcNotImplemented()); + template + void - // first build the matrices M and Q - // described in the documentation - FullMatrix M (fe.dofs_per_cell, fe.dofs_per_cell); - FullMatrix Q (fe.dofs_per_cell, rhs_quadrature.size()); + compute_projection_from_quadrature_points_matrix (const FiniteElement &fe, + const Quadrature &lhs_quadrature, + const Quadrature &rhs_quadrature, + FullMatrix &X) + { + Assert (fe.n_components() == 1, ExcNotImplemented()); + + // first build the matrices M and Q + // described in the documentation + FullMatrix M (fe.dofs_per_cell, fe.dofs_per_cell); + FullMatrix Q (fe.dofs_per_cell, rhs_quadrature.size()); + + for (unsigned int i=0; i M_inverse (fe.dofs_per_cell, fe.dofs_per_cell); + M_inverse.invert (M); + + // finally compute the result + X.reinit (fe.dofs_per_cell, rhs_quadrature.size()); + M_inverse.mmult (X, Q); + + Assert (X.m() == fe.dofs_per_cell, ExcInternalError()); + Assert (X.n() == rhs_quadrature.size(), ExcInternalError()); + } - for (unsigned int i=0; i + void + compute_interpolation_to_quadrature_points_matrix (const FiniteElement &fe, + const Quadrature &quadrature, + FullMatrix &I_q) + { + Assert (fe.n_components() == 1, ExcNotImplemented()); + Assert (I_q.m() == quadrature.size(), + ExcMessage ("Wrong matrix size")); + Assert (I_q.n() == fe.dofs_per_cell, ExcMessage ("Wrong matrix size")); + + for (unsigned int q=0; q M_inverse (fe.dofs_per_cell, fe.dofs_per_cell); - M_inverse.invert (M); - // finally compute the result - X.reinit (fe.dofs_per_cell, rhs_quadrature.size()); - M_inverse.mmult (X, Q); - Assert (X.m() == fe.dofs_per_cell, ExcInternalError()); - Assert (X.n() == rhs_quadrature.size(), ExcInternalError()); - } + template + void + compute_projection_from_quadrature_points( + const FullMatrix &projection_matrix, + const std::vector< Tensor<1, dim > > &vector_of_tensors_at_qp, + std::vector< Tensor<1, dim > > &vector_of_tensors_at_nodes) + { + // check that the number columns of the projection_matrix + // matches the size of the vector_of_tensors_at_qp + Assert(projection_matrix.n_cols() == vector_of_tensors_at_qp.size(), + ExcDimensionMismatch(projection_matrix.n_cols(), + vector_of_tensors_at_qp.size())); + + // check that the number rows of the projection_matrix + // matches the size of the vector_of_tensors_at_nodes + Assert(projection_matrix.n_rows() == vector_of_tensors_at_nodes.size(), + ExcDimensionMismatch(projection_matrix.n_rows(), + vector_of_tensors_at_nodes.size())); + + // number of support points (nodes) to project to + const unsigned int n_support_points = projection_matrix.n_rows(); + // number of quadrature points to project from + const unsigned int n_quad_points = projection_matrix.n_cols(); + + // component projected to the nodes + Vector component_at_node(n_support_points); + // component at the quadrature point + Vector component_at_qp(n_quad_points); + + for (unsigned int ii = 0; ii < dim; ++ii) + { + component_at_qp = 0; - template - void - compute_interpolation_to_quadrature_points_matrix (const FiniteElement &fe, - const Quadrature &quadrature, - FullMatrix &I_q) - { - Assert (fe.n_components() == 1, ExcNotImplemented()); - Assert (I_q.m() == quadrature.size(), - ExcMessage ("Wrong matrix size")); - Assert (I_q.n() == fe.dofs_per_cell, ExcMessage ("Wrong matrix size")); - - for (unsigned int q=0; q nodes + // component_at_node = projection_matrix_u * component_at_qp + projection_matrix.vmult(component_at_node, component_at_qp); + // rewrite the projection of the components + // back into the vector of tensors + for (unsigned int nn =0; nn - void - compute_projection_from_quadrature_points( - const FullMatrix &projection_matrix, - const std::vector< Tensor<1, dim > > &vector_of_tensors_at_qp, - std::vector< Tensor<1, dim > > &vector_of_tensors_at_nodes) - { - // check that the number columns of the projection_matrix - // matches the size of the vector_of_tensors_at_qp - Assert(projection_matrix.n_cols() == vector_of_tensors_at_qp.size(), - ExcDimensionMismatch(projection_matrix.n_cols(), - vector_of_tensors_at_qp.size())); - - // check that the number rows of the projection_matrix - // matches the size of the vector_of_tensors_at_nodes - Assert(projection_matrix.n_rows() == vector_of_tensors_at_nodes.size(), - ExcDimensionMismatch(projection_matrix.n_rows(), - vector_of_tensors_at_nodes.size())); - - // number of support points (nodes) to project to - const unsigned int n_support_points = projection_matrix.n_rows(); - // number of quadrature points to project from - const unsigned int n_quad_points = projection_matrix.n_cols(); - - // component projected to the nodes - Vector component_at_node(n_support_points); - // component at the quadrature point - Vector component_at_qp(n_quad_points); - - for (unsigned int ii = 0; ii < dim; ++ii) - { - component_at_qp = 0; + template + void + compute_projection_from_quadrature_points( + const FullMatrix &projection_matrix, + const std::vector< SymmetricTensor<2, dim > > &vector_of_tensors_at_qp, + std::vector< SymmetricTensor<2, dim > > &vector_of_tensors_at_nodes) + { - // populate the vector of components at the qps - // from vector_of_tensors_at_qp - // vector_of_tensors_at_qp data is in form: - // columns: 0, 1, ..., dim - // rows: 0,1,...., n_quad_points - // so extract the ii'th column of vector_of_tensors_at_qp - for (unsigned int q = 0; q < n_quad_points; ++q) - { - component_at_qp(q) = vector_of_tensors_at_qp[q][ii]; - } + // check that the number columns of the projection_matrix + // matches the size of the vector_of_tensors_at_qp + Assert(projection_matrix.n_cols() == vector_of_tensors_at_qp.size(), + ExcDimensionMismatch(projection_matrix.n_cols(), + vector_of_tensors_at_qp.size())); + + // check that the number rows of the projection_matrix + // matches the size of the vector_of_tensors_at_nodes + Assert(projection_matrix.n_rows() == vector_of_tensors_at_nodes.size(), + ExcDimensionMismatch(projection_matrix.n_rows(), + vector_of_tensors_at_nodes.size())); + + // number of support points (nodes) + const unsigned int n_support_points = projection_matrix.n_rows(); + // number of quadrature points to project from + const unsigned int n_quad_points = projection_matrix.n_cols(); + + // number of unique entries in a symmetric second-order tensor + const unsigned int n_independent_components = + SymmetricTensor<2, dim >::n_independent_components; + + // component projected to the nodes + Vector component_at_node(n_support_points); + // component at the quadrature point + Vector component_at_qp(n_quad_points); + + // loop over the number of unique dimensions of the tensor + for (unsigned int ii = 0; ii < n_independent_components; ++ii) + { - // project from the qps -> nodes - // component_at_node = projection_matrix_u * component_at_qp - projection_matrix.vmult(component_at_node, component_at_qp); + component_at_qp = 0; - // rewrite the projection of the components - // back into the vector of tensors - for (unsigned int nn =0; nn row_column_index = SymmetricTensor< 2, dim >::unrolled_to_component_indices(ii); + const unsigned int row = row_column_index[0]; + const unsigned int column = row_column_index[1]; + // populate the vector of components at the qps + // from vector_of_tensors_at_qp + // vector_of_tensors_at_qp is in form: + // columns: 0, 1, ..., n_independent_components + // rows: 0,1,...., n_quad_points + // so extract the ii'th column of vector_of_tensors_at_qp + for (unsigned int q = 0; q < n_quad_points; ++q) + { + component_at_qp(q) = (vector_of_tensors_at_qp[q])[row][column]; + } + // project from the qps -> nodes + // component_at_node = projection_matrix_u * component_at_qp + projection_matrix.vmult(component_at_node, component_at_qp); - template - void - compute_projection_from_quadrature_points( - const FullMatrix &projection_matrix, - const std::vector< SymmetricTensor<2, dim > > &vector_of_tensors_at_qp, - std::vector< SymmetricTensor<2, dim > > &vector_of_tensors_at_nodes) - { + // rewrite the projection of the components back into the vector of tensors + for (unsigned int nn =0; nn ::n_independent_components; - - // component projected to the nodes - Vector component_at_node(n_support_points); - // component at the quadrature point - Vector component_at_qp(n_quad_points); - - // loop over the number of unique dimensions of the tensor - for (unsigned int ii = 0; ii < n_independent_components; ++ii) - { - component_at_qp = 0; - // row-column entry of tensor corresponding the unrolled index - TableIndices<2> row_column_index = SymmetricTensor< 2, dim >::unrolled_to_component_indices(ii); - const unsigned int row = row_column_index[0]; - const unsigned int column = row_column_index[1]; + template + void + compute_projection_from_face_quadrature_points_matrix (const FiniteElement &fe, + const Quadrature &lhs_quadrature, + const Quadrature &rhs_quadrature, + const typename DoFHandler::active_cell_iterator &cell, + const unsigned int face, + FullMatrix &X) + { + Assert (fe.n_components() == 1, ExcNotImplemented()); + Assert (lhs_quadrature.size () > fe.degree, ExcNotGreaterThan (lhs_quadrature.size (), fe.degree)); - // populate the vector of components at the qps - // from vector_of_tensors_at_qp - // vector_of_tensors_at_qp is in form: - // columns: 0, 1, ..., n_independent_components - // rows: 0,1,...., n_quad_points - // so extract the ii'th column of vector_of_tensors_at_qp - for (unsigned int q = 0; q < n_quad_points; ++q) - { - component_at_qp(q) = (vector_of_tensors_at_qp[q])[row][column]; - } - // project from the qps -> nodes - // component_at_node = projection_matrix_u * component_at_qp - projection_matrix.vmult(component_at_node, component_at_qp); - // rewrite the projection of the components back into the vector of tensors - for (unsigned int nn =0; nn M (fe.dofs_per_cell, fe.dofs_per_cell); + FullMatrix Q (fe.dofs_per_cell, rhs_quadrature.size()); + + { + // need an FEFaceValues object to evaluate shape function + // values on the specified face. + FEFaceValues fe_face_values (fe, lhs_quadrature, update_values); + fe_face_values.reinit (cell, face); // setup shape_value on this face. + + for (unsigned int i=0; i fe_face_values (fe, rhs_quadrature, update_values); + fe_face_values.reinit (cell, face); // setup shape_value on this face. + for (unsigned int i=0; i M_inverse (fe.dofs_per_cell, fe.dofs_per_cell); + M_inverse.invert (M); - template - void - compute_projection_from_face_quadrature_points_matrix (const FiniteElement &fe, - const Quadrature &lhs_quadrature, - const Quadrature &rhs_quadrature, - const typename DoFHandler::active_cell_iterator &cell, - const unsigned int face, - FullMatrix &X) - { - Assert (fe.n_components() == 1, ExcNotImplemented()); - Assert (lhs_quadrature.size () > fe.degree, ExcNotGreaterThan (lhs_quadrature.size (), fe.degree)); + // finally compute the result + X.reinit (fe.dofs_per_cell, rhs_quadrature.size()); + M_inverse.mmult (X, Q); + Assert (X.m() == fe.dofs_per_cell, ExcInternalError()); + Assert (X.n() == rhs_quadrature.size(), ExcInternalError()); + } - // build the matrices M and Q - // described in the documentation - FullMatrix M (fe.dofs_per_cell, fe.dofs_per_cell); - FullMatrix Q (fe.dofs_per_cell, rhs_quadrature.size()); + template + void + hierarchic_to_lexicographic_numbering (unsigned int degree, std::vector &h2l) + { + // number of support points in each + // direction + const unsigned int n = degree+1; + + unsigned int dofs_per_cell = n; + for (unsigned int i=1; i fe_face_values (fe, lhs_quadrature, update_values); - fe_face_values.reinit (cell, face); // setup shape_value on this face. - - for (unsigned int i=0; i fe_face_values (fe, rhs_quadrature, update_values); - fe_face_values.reinit (cell, face); // setup shape_value on this face. + h2l[0] = 0; + h2l[1] = dofs_per_cell-1; + for (unsigned int i=2; i M_inverse (fe.dofs_per_cell, fe.dofs_per_cell); - M_inverse.invert (M); - - // finally compute the result - X.reinit (fe.dofs_per_cell, rhs_quadrature.size()); - M_inverse.mmult (X, Q); - - Assert (X.m() == fe.dofs_per_cell, ExcInternalError()); - Assert (X.n() == rhs_quadrature.size(), ExcInternalError()); - } - - - - template - void - hierarchic_to_lexicographic_numbering (unsigned int degree, std::vector &h2l) - { - // number of support points in each - // direction - const unsigned int n = degree+1; - - unsigned int dofs_per_cell = n; - for (unsigned int i=1; i - void - hierarchic_to_lexicographic_numbering (const FiniteElementData &fe, - std::vector &h2l) - { - Assert (h2l.size() == fe.dofs_per_cell, - ExcDimensionMismatch (h2l.size(), fe.dofs_per_cell)); - hierarchic_to_lexicographic_numbering (fe.dofs_per_line+1, h2l); - } + template + void + hierarchic_to_lexicographic_numbering (const FiniteElementData &fe, + std::vector &h2l) + { + Assert (h2l.size() == fe.dofs_per_cell, + ExcDimensionMismatch (h2l.size(), fe.dofs_per_cell)); + hierarchic_to_lexicographic_numbering (fe.dofs_per_line+1, h2l); + } - template - std::vector - hierarchic_to_lexicographic_numbering (const FiniteElementData &fe) - { - Assert (fe.n_components() == 1, ExcInvalidFE()); - std::vector h2l(fe.dofs_per_cell); - hierarchic_to_lexicographic_numbering (fe.dofs_per_line+1, h2l); - return (h2l); - } + template + std::vector + hierarchic_to_lexicographic_numbering (const FiniteElementData &fe) + { + Assert (fe.n_components() == 1, ExcInvalidFE()); + std::vector h2l(fe.dofs_per_cell); + hierarchic_to_lexicographic_numbering (fe.dofs_per_line+1, h2l); + return (h2l); + } - template - void - lexicographic_to_hierarchic_numbering (const FiniteElementData &fe, - std::vector &l2h) - { - l2h = lexicographic_to_hierarchic_numbering (fe); - } + template + void + lexicographic_to_hierarchic_numbering (const FiniteElementData &fe, + std::vector &l2h) + { + l2h = lexicographic_to_hierarchic_numbering (fe); + } - template - std::vector - lexicographic_to_hierarchic_numbering (const FiniteElementData &fe) - { - return Utilities::invert_permutation(hierarchic_to_lexicographic_numbering (fe)); - } + template + std::vector + lexicographic_to_hierarchic_numbering (const FiniteElementData &fe) + { + return Utilities::invert_permutation(hierarchic_to_lexicographic_numbering (fe)); + } - } // end of namespace FETools +} // end of namespace FETools - /*-------------- Explicit Instantiations -------------------------------*/ +/*-------------- Explicit Instantiations -------------------------------*/ #include "fe_tools.inst" - /*---------------------------- fe_tools.cc ---------------------------*/ +/*---------------------------- fe_tools.cc ---------------------------*/ - DEAL_II_NAMESPACE_CLOSE +DEAL_II_NAMESPACE_CLOSE -- 2.39.5