From: Luca Heltai Date: Tue, 1 Mar 2016 18:22:44 +0000 (+0100) Subject: Removed some old stuff. X-Git-Tag: v8.5.0-rc1~1131^2~34 X-Git-Url: https://gitweb.dealii.org/cgi-bin/gitweb.cgi?a=commitdiff_plain;h=21f037b4c0bf679f72161899165f3b3f39ecb5cb;p=dealii.git Removed some old stuff. --- diff --git a/include/deal.II/fe/mapping_manifold.h b/include/deal.II/fe/mapping_manifold.h index 7819542107..31beb2326e 100644 --- a/include/deal.II/fe/mapping_manifold.h +++ b/include/deal.II/fe/mapping_manifold.h @@ -196,84 +196,15 @@ public: const Quadrature &quadrature, const unsigned int n_original_q_points); - /** - * Compute the values and/or derivatives of the shape functions used for - * the mapping. - * - * Which values, derivatives, or higher order derivatives are computed is - * determined by which of the member arrays have nonzero sizes. They are - * typically set to their appropriate sizes by the initialize() and - * initialize_face() functions, which indeed call this function - * internally. However, it is possible (and at times useful) to do the - * resizing by hand and then call this function directly. An example is in - * a Newton iteration where we update the location of a quadrature point - * (e.g., in MappingQ::transform_real_to_uni_cell()) and need to re- - * compute the mapping and its derivatives at this location, but have - * already sized all internal arrays correctly. - */ - void compute_shape_function_values (const std::vector > &unit_points); - - - /** - * Shape function at quadrature point. Shape functions are in tensor - * product order, so vertices must be reordered to obtain transformation. - */ - const double &shape (const unsigned int qpoint, - const unsigned int shape_nr) const; - - /** - * Shape function at quadrature point. See above. - */ - double &shape (const unsigned int qpoint, - const unsigned int shape_nr); - - /** - * Gradient of shape function in quadrature point. See above. - */ - const Tensor<1,dim> &derivative (const unsigned int qpoint, - const unsigned int shape_nr) const; - - /** - * Gradient of shape function in quadrature point. See above. - */ - Tensor<1,dim> &derivative (const unsigned int qpoint, - const unsigned int shape_nr); - - /** - * Second derivative of shape function in quadrature point. See above. - */ - const Tensor<2,dim> &second_derivative (const unsigned int qpoint, - const unsigned int shape_nr) const; - - /** - * Second derivative of shape function in quadrature point. See above. - */ - Tensor<2,dim> &second_derivative (const unsigned int qpoint, - const unsigned int shape_nr); /** - * third derivative of shape function in quadrature point. See above. + * Compute the weights associated to the Manifold object, that + * need to be passed when computing the location of the quadrature + * points. */ - const Tensor<3,dim> &third_derivative (const unsigned int qpoint, - const unsigned int shape_nr) const; - - /** - * third derivative of shape function in quadrature point. See above. - */ - Tensor<3,dim> &third_derivative (const unsigned int qpoint, - const unsigned int shape_nr); - - /** - * fourth derivative of shape function in quadrature point. See above. - */ - const Tensor<4,dim> &fourth_derivative (const unsigned int qpoint, - const unsigned int shape_nr) const; + void + compute_manifold_quadrature_weights(const Quadrature &quadrature); - /** - * fourth derivative of shape function in quadrature point. See above. - */ - Tensor<4,dim> &fourth_derivative (const unsigned int qpoint, - const unsigned int shape_nr); /** * Return an estimate (in bytes) or the memory consumption of this object. @@ -281,116 +212,123 @@ public: virtual std::size_t memory_consumption () const; /** - * Values of shape functions. Access by function @p shape. - * - * Computed once. - */ - std::vector shape_values; - - /** - * Values of shape function derivatives. Access by function @p derivative. + * Values of manifold quadrature formulas. * - * Computed once. + * Computed each. */ - std::vector > shape_derivatives; + std::vector > cell_manifold_quadratures; /** - * Values of shape function second derivatives. Access by function @p - * second_derivative. + * Values of quadrature weights for manifold quadrature formulas. * * Computed once. */ - std::vector > shape_second_derivatives; - - /** - * Values of shape function third derivatives. Access by function @p - * second_derivative. - * - * Computed once. - */ - std::vector > shape_third_derivatives; - - /** - * Values of shape function fourth derivatives. Access by function @p - * second_derivative. - * - * Computed once. - */ - std::vector > shape_fourth_derivatives; - - /** - * Unit tangential vectors. Used for the computation of boundary forms and - * normal vectors. - * - * This vector has (dim-1)GeometryInfo::faces_per_cell entries. The first - * GeometryInfo::faces_per_cell contain the vectors in the first - * tangential direction for each face; the second set of - * GeometryInfo::faces_per_cell entries contain the vectors in the second - * tangential direction (only in 3d, since there we have 2 tangential - * directions per face), etc. - * - * Filled once. - */ - std::vector > > unit_tangentials; - - /** - * The polynomial degree of the mapping. Since the objects here are also - * used (with minor adjustments) by MappingQ, we need to store this. - */ - unsigned int polynomial_degree; - - /** - * Number of shape functions. If this is a Q1 mapping, then it is simply - * the number of vertices per cell. However, since also derived classes - * use this class (e.g. the Mapping_Q() class), the number of shape - * functions may also be different. - * - * In general, it is $(p+1)^\text{dim}$, where $p$ is the polynomial - * degree of the mapping. - */ - const unsigned int n_shape_functions; - - /** - * Tensors of covariant transformation at each of the quadrature points. - * The matrix stored is the Jacobian * G^{-1}, where G = Jacobian^{t} * - * Jacobian, is the first fundamental form of the map; if dim=spacedim - * then it reduces to the transpose of the inverse of the Jacobian matrix, - * which itself is stored in the @p contravariant field of this structure. - * - * Computed on each cell. - */ - mutable std::vector > covariant; - - /** - * Tensors of contravariant transformation at each of the quadrature - * points. The contravariant matrix is the Jacobian of the transformation, - * i.e. $J_{ij}=dx_i/d\hat x_j$. - * - * Computed on each cell. - */ - mutable std::vector< DerivativeForm<1,dim,spacedim> > contravariant; - - /** - * Auxiliary vectors for internal use. - */ - mutable std::vector > > aux; - - /** - * Stores the support points of the mapping shape functions on the @p - * cell_of_current_support_points. - */ - mutable std::vector > mapping_support_points; - - /** - * Stores the cell of which the @p mapping_support_points are stored. - */ - mutable typename Triangulation::cell_iterator cell_of_current_support_points; - - /** - * The determinant of the Jacobian in each quadrature point. Filled if - * #update_volume_elements. - */ - mutable std::vector volume_elements; + std::vector > cell_manifold_quadratures_weights; + + // /** + // * Values of shape function derivatives. Access by function @p derivative. + // * + // * Computed once. + // */ + // std::vector > shape_derivatives; + + // /** + // * Values of shape function second derivatives. Access by function @p + // * second_derivative. + // * + // * Computed once. + // */ + // std::vector > shape_second_derivatives; + + // /** + // * Values of shape function third derivatives. Access by function @p + // * second_derivative. + // * + // * Computed once. + // */ + // std::vector > shape_third_derivatives; + + // /** + // * Values of shape function fourth derivatives. Access by function @p + // * second_derivative. + // * + // * Computed once. + // */ + // std::vector > shape_fourth_derivatives; + + // /** + // * Unit tangential vectors. Used for the computation of boundary forms and + // * normal vectors. + // * + // * This vector has (dim-1)GeometryInfo::faces_per_cell entries. The first + // * GeometryInfo::faces_per_cell contain the vectors in the first + // * tangential direction for each face; the second set of + // * GeometryInfo::faces_per_cell entries contain the vectors in the second + // * tangential direction (only in 3d, since there we have 2 tangential + // * directions per face), etc. + // * + // * Filled once. + // */ + // std::vector > > unit_tangentials; + + // /** + // * The polynomial degree of the mapping. Since the objects here are also + // * used (with minor adjustments) by MappingQ, we need to store this. + // */ + // unsigned int polynomial_degree; + + // /** + // * Number of shape functions. If this is a Q1 mapping, then it is simply + // * the number of vertices per cell. However, since also derived classes + // * use this class (e.g. the Mapping_Q() class), the number of shape + // * functions may also be different. + // * + // * In general, it is $(p+1)^\text{dim}$, where $p$ is the polynomial + // * degree of the mapping. + // */ + // const unsigned int n_shape_functions; + + // /** + // * Tensors of covariant transformation at each of the quadrature points. + // * The matrix stored is the Jacobian * G^{-1}, where G = Jacobian^{t} * + // * Jacobian, is the first fundamental form of the map; if dim=spacedim + // * then it reduces to the transpose of the inverse of the Jacobian matrix, + // * which itself is stored in the @p contravariant field of this structure. + // * + // * Computed on each cell. + // */ + // mutable std::vector > covariant; + + // /** + // * Tensors of contravariant transformation at each of the quadrature + // * points. The contravariant matrix is the Jacobian of the transformation, + // * i.e. $J_{ij}=dx_i/d\hat x_j$. + // * + // * Computed on each cell. + // */ + // mutable std::vector< DerivativeForm<1,dim,spacedim> > contravariant; + + // /** + // * Auxiliary vectors for internal use. + // */ + // mutable std::vector > > aux; + + // /** + // * Stores the support points of the mapping shape functions on the @p + // * cell_of_current_support_points. + // */ + // mutable std::vector > mapping_support_points; + + // /** + // * Stores the cell of which the @p mapping_support_points are stored. + // */ + // mutable typename Triangulation::cell_iterator cell_of_current_support_points; + + // /** + // * The determinant of the Jacobian in each quadrature point. Filled if + // * #update_volume_elements. + // */ + // mutable std::vector volume_elements; }; @@ -450,28 +388,9 @@ public: protected: /** - * The degree of the polynomials used as shape functions for the mapping of - * cells. - */ - const unsigned int polynomial_degree; - - /* - * The default line support points. These are used when computing - * the location in real space of the support points on lines and - * quads, which are asked to the Manifold class. - * - * The number of quadrature points depends on the degree of this - * class, and it matches the number of degrees of freedom of an - * FE_Q<1>(this->degree). - */ - QGaussLobatto<1> line_support_points; - - /** - * An FE_Q object which is only needed in 3D, since it knows how to reorder - * shape functions/DoFs on non-standard faces. This is used to reorder - * support points in the same way. + * An FE_Q object, used to compute weights for Manifold quadratures. */ - const std_cxx11::unique_ptr > fe_q; + const FE_Q fe_q; /** * A table of weights by which we multiply the locations of the support @@ -582,6 +501,7 @@ protected: * functions on its MappingManifold(1) sub-object. */ template friend class MappingQ; + }; @@ -592,134 +512,124 @@ protected: #ifndef DOXYGEN -template -inline -const double & -MappingManifold::InternalData::shape (const unsigned int qpoint, - const unsigned int shape_nr) const -{ - Assert(qpoint*n_shape_functions + shape_nr < shape_values.size(), - ExcIndexRange(qpoint*n_shape_functions + shape_nr, 0, - shape_values.size())); - return shape_values [qpoint*n_shape_functions + shape_nr]; -} - - - -template -inline -double & -MappingManifold::InternalData::shape (const unsigned int qpoint, - const unsigned int shape_nr) -{ - Assert(qpoint*n_shape_functions + shape_nr < shape_values.size(), - ExcIndexRange(qpoint*n_shape_functions + shape_nr, 0, - shape_values.size())); - return shape_values [qpoint*n_shape_functions + shape_nr]; -} - - -template -inline -const Tensor<1,dim> & -MappingManifold::InternalData::derivative (const unsigned int qpoint, - const unsigned int shape_nr) const -{ - Assert(qpoint*n_shape_functions + shape_nr < shape_derivatives.size(), - ExcIndexRange(qpoint*n_shape_functions + shape_nr, 0, - shape_derivatives.size())); - return shape_derivatives [qpoint*n_shape_functions + shape_nr]; -} - - - -template -inline -Tensor<1,dim> & -MappingManifold::InternalData::derivative (const unsigned int qpoint, - const unsigned int shape_nr) -{ - Assert(qpoint*n_shape_functions + shape_nr < shape_derivatives.size(), - ExcIndexRange(qpoint*n_shape_functions + shape_nr, 0, - shape_derivatives.size())); - return shape_derivatives [qpoint*n_shape_functions + shape_nr]; -} - - -template -inline -const Tensor<2,dim> & -MappingManifold::InternalData::second_derivative (const unsigned int qpoint, - const unsigned int shape_nr) const -{ - Assert(qpoint*n_shape_functions + shape_nr < shape_second_derivatives.size(), - ExcIndexRange(qpoint*n_shape_functions + shape_nr, 0, - shape_second_derivatives.size())); - return shape_second_derivatives [qpoint*n_shape_functions + shape_nr]; -} - - -template -inline -Tensor<2,dim> & -MappingManifold::InternalData::second_derivative (const unsigned int qpoint, - const unsigned int shape_nr) -{ - Assert(qpoint*n_shape_functions + shape_nr < shape_second_derivatives.size(), - ExcIndexRange(qpoint*n_shape_functions + shape_nr, 0, - shape_second_derivatives.size())); - return shape_second_derivatives [qpoint*n_shape_functions + shape_nr]; -} - -template -inline -const Tensor<3,dim> & -MappingManifold::InternalData::third_derivative (const unsigned int qpoint, - const unsigned int shape_nr) const -{ - Assert(qpoint*n_shape_functions + shape_nr < shape_third_derivatives.size(), - ExcIndexRange(qpoint*n_shape_functions + shape_nr, 0, - shape_third_derivatives.size())); - return shape_third_derivatives [qpoint*n_shape_functions + shape_nr]; -} - - -template -inline -Tensor<3,dim> & -MappingManifold::InternalData::third_derivative (const unsigned int qpoint, - const unsigned int shape_nr) -{ - Assert(qpoint*n_shape_functions + shape_nr < shape_third_derivatives.size(), - ExcIndexRange(qpoint*n_shape_functions + shape_nr, 0, - shape_third_derivatives.size())); - return shape_third_derivatives [qpoint*n_shape_functions + shape_nr]; -} - - -template -inline -const Tensor<4,dim> & -MappingManifold::InternalData::fourth_derivative (const unsigned int qpoint, - const unsigned int shape_nr) const -{ - Assert(qpoint*n_shape_functions + shape_nr < shape_fourth_derivatives.size(), - ExcIndexRange(qpoint*n_shape_functions + shape_nr, 0, - shape_fourth_derivatives.size())); - return shape_fourth_derivatives [qpoint*n_shape_functions + shape_nr]; -} +// template +// inline +// const Tensor<1,dim> & +// MappingManifold::InternalData::derivative (const unsigned int qpoint, +// const unsigned int shape_nr) const +// { +// Assert(qpoint*n_shape_functions + shape_nr < shape_derivatives.size(), +// ExcIndexRange(qpoint*n_shape_functions + shape_nr, 0, +// shape_derivatives.size())); +// return shape_derivatives [qpoint*n_shape_functions + shape_nr]; +// } + + + +// template +// inline +// Tensor<1,dim> & +// MappingManifold::InternalData::derivative (const unsigned int qpoint, +// const unsigned int shape_nr) +// { +// Assert(qpoint*n_shape_functions + shape_nr < shape_derivatives.size(), +// ExcIndexRange(qpoint*n_shape_functions + shape_nr, 0, +// shape_derivatives.size())); +// return shape_derivatives [qpoint*n_shape_functions + shape_nr]; +// } + + +// template +// inline +// const Tensor<2,dim> & +// MappingManifold::InternalData::second_derivative (const unsigned int qpoint, +// const unsigned int shape_nr) const +// { +// Assert(qpoint*n_shape_functions + shape_nr < shape_second_derivatives.size(), +// ExcIndexRange(qpoint*n_shape_functions + shape_nr, 0, +// shape_second_derivatives.size())); +// return shape_second_derivatives [qpoint*n_shape_functions + shape_nr]; +// } + + +// template +// inline +// Tensor<2,dim> & +// MappingManifold::InternalData::second_derivative (const unsigned int qpoint, +// const unsigned int shape_nr) +// { +// Assert(qpoint*n_shape_functions + shape_nr < shape_second_derivatives.size(), +// ExcIndexRange(qpoint*n_shape_functions + shape_nr, 0, +// shape_second_derivatives.size())); +// return shape_second_derivatives [qpoint*n_shape_functions + shape_nr]; +// } + +// template +// inline +// const Tensor<3,dim> & +// MappingManifold::InternalData::third_derivative (const unsigned int qpoint, +// const unsigned int shape_nr) const +// { +// Assert(qpoint*n_shape_functions + shape_nr < shape_third_derivatives.size(), +// ExcIndexRange(qpoint*n_shape_functions + shape_nr, 0, +// shape_third_derivatives.size())); +// return shape_third_derivatives [qpoint*n_shape_functions + shape_nr]; +// } + + +// template +// inline +// Tensor<3,dim> & +// MappingManifold::InternalData::third_derivative (const unsigned int qpoint, +// const unsigned int shape_nr) +// { +// Assert(qpoint*n_shape_functions + shape_nr < shape_third_derivatives.size(), +// ExcIndexRange(qpoint*n_shape_functions + shape_nr, 0, +// shape_third_derivatives.size())); +// return shape_third_derivatives [qpoint*n_shape_functions + shape_nr]; +// } + + +// template +// inline +// const Tensor<4,dim> & +// MappingManifold::InternalData::fourth_derivative (const unsigned int qpoint, +// const unsigned int shape_nr) const +// { +// Assert(qpoint*n_shape_functions + shape_nr < shape_fourth_derivatives.size(), +// ExcIndexRange(qpoint*n_shape_functions + shape_nr, 0, +// shape_fourth_derivatives.size())); +// return shape_fourth_derivatives [qpoint*n_shape_functions + shape_nr]; +// } + + +// template +// inline +// Tensor<4,dim> & +// MappingManifold::InternalData::fourth_derivative (const unsigned int qpoint, +// const unsigned int shape_nr) +// { +// Assert(qpoint*n_shape_functions + shape_nr < shape_fourth_derivatives.size(), +// ExcIndexRange(qpoint*n_shape_functions + shape_nr, 0, +// shape_fourth_derivatives.size())); +// return shape_fourth_derivatives [qpoint*n_shape_functions + shape_nr]; +// } template inline -Tensor<4,dim> & -MappingManifold::InternalData::fourth_derivative (const unsigned int qpoint, - const unsigned int shape_nr) +void +MappingManifold::InternalData::compute_manifold_quadrature_weights (const Quadrature &quad) { - Assert(qpoint*n_shape_functions + shape_nr < shape_fourth_derivatives.size(), - ExcIndexRange(qpoint*n_shape_functions + shape_nr, 0, - shape_fourth_derivatives.size())); - return shape_fourth_derivatives [qpoint*n_shape_functions + shape_nr]; + static FE_Q fe_q(1); + cell_manifold_quadratures_weights.resize(quad.size(), std::vector(GeometryInfo::vertices_per_cell)); + for (unsigned int q=0; q::vertices_per_cell; ++i) + { + cell_manifold_quadratures_weights[q][i] = fe_q.shape_value(i, quad.point(q)); + } + } } diff --git a/source/fe/mapping_manifold.cc b/source/fe/mapping_manifold.cc index 5aa00c13a4..78796b8ff0 100644 --- a/source/fe/mapping_manifold.cc +++ b/source/fe/mapping_manifold.cc @@ -40,632 +40,8 @@ DEAL_II_NAMESPACE_OPEN -namespace internal -{ - namespace MappingQ1 - { - namespace - { - - // These are left as templates on the spatial dimension (even though dim - // == spacedim must be true for them to make sense) because templates are - // expanded before the compiler eliminates code due to the 'if (dim == - // spacedim)' statement (see the body of the general - // transform_real_to_unit_cell). - template - Point<1> - transform_real_to_unit_cell - (const std_cxx11::array, GeometryInfo<1>::vertices_per_cell> &vertices, - const Point &p) - { - Assert(spacedim == 1, ExcInternalError()); - return Point<1>((p[0] - vertices[0](0))/(vertices[1](0) - vertices[0](0))); - } - - - - template - Point<2> - transform_real_to_unit_cell - (const std_cxx11::array, GeometryInfo<2>::vertices_per_cell> &vertices, - const Point &p) - { - Assert(spacedim == 2, ExcInternalError()); - const double x = p(0); - const double y = p(1); - - const double x0 = vertices[0](0); - const double x1 = vertices[1](0); - const double x2 = vertices[2](0); - const double x3 = vertices[3](0); - - const double y0 = vertices[0](1); - const double y1 = vertices[1](1); - const double y2 = vertices[2](1); - const double y3 = vertices[3](1); - - const double a = (x1 - x3)*(y0 - y2) - (x0 - x2)*(y1 - y3); - const double b = -(x0 - x1 - x2 + x3)*y + (x - 2*x1 + x3)*y0 - (x - 2*x0 + x2)*y1 - - (x - x1)*y2 + (x - x0)*y3; - const double c = (x0 - x1)*y - (x - x1)*y0 + (x - x0)*y1; - - const double discriminant = b*b - 4*a*c; - // exit if the point is not in the cell (this is the only case where the - // discriminant is negative) - if (discriminant < 0.0) - { - AssertThrow (false, - (typename Mapping::ExcTransformationFailed())); - } - - double eta1; - double eta2; - // special case #1: if a is zero, then use the linear formula - if (a == 0.0 && b != 0.0) - { - eta1 = -c/b; - eta2 = -c/b; - } - // special case #2: if c is very small or the square root of the - // discriminant is nearly b. - else if (std::abs(c) < 1e-12*std::abs(b) - || std::abs(std::sqrt(discriminant) - b) <= 1e-14*std::abs(b)) - { - eta1 = (-b - std::sqrt(discriminant)) / (2*a); - eta2 = (-b + std::sqrt(discriminant)) / (2*a); - } - // finally, use the numerically stable version of the quadratic formula: - else - { - eta1 = 2*c / (-b - std::sqrt(discriminant)); - eta2 = 2*c / (-b + std::sqrt(discriminant)); - } - // pick the one closer to the center of the cell. - const double eta = (std::abs(eta1 - 0.5) < std::abs(eta2 - 0.5)) ? eta1 : eta2; - - /* - * There are two ways to compute xi from eta, but either one may have a - * zero denominator. - */ - const double subexpr0 = -eta*x2 + x0*(eta - 1); - const double xi_denominator0 = eta*x3 - x1*(eta - 1) + subexpr0; - const double max_x = std::max(std::max(std::abs(x0), std::abs(x1)), - std::max(std::abs(x2), std::abs(x3))); - - if (std::abs(xi_denominator0) > 1e-10*max_x) - { - const double xi = (x + subexpr0)/xi_denominator0; - return Point<2>(xi, eta); - } - else - { - const double max_y = std::max(std::max(std::abs(y0), std::abs(y1)), - std::max(std::abs(y2), std::abs(y3))); - const double subexpr1 = -eta*y2 + y0*(eta - 1); - const double xi_denominator1 = eta*y3 - y1*(eta - 1) + subexpr1; - if (std::abs(xi_denominator1) > 1e-10*max_y) - { - const double xi = (subexpr1 + y)/xi_denominator1; - return Point<2>(xi, eta); - } - else // give up and try Newton iteration - { - AssertThrow (false, - (typename Mapping::ExcTransformationFailed())); - } - } - // bogus return to placate compiler. It should not be possible to get - // here. - Assert(false, ExcInternalError()); - return Point<2>(std::numeric_limits::quiet_NaN(), - std::numeric_limits::quiet_NaN()); - } - - - - template - Point<3> - transform_real_to_unit_cell - (const std_cxx11::array, GeometryInfo<3>::vertices_per_cell> &/*vertices*/, - const Point &/*p*/) - { - // It should not be possible to get here - Assert(false, ExcInternalError()); - return Point<3>(); - } - - - - /** - * Compute an initial guess to pass to the Newton method in - * transform_real_to_unit_cell. For the initial guess we proceed in the - * following way: - *
    - *
  • find the least square dim-dimensional plane approximating the cell - * vertices, i.e. we find an affine map A x_hat + b from the reference cell - * to the real space. - *
  • Solve the equation A x_hat + b = p for x_hat - *
  • This x_hat is the initial solution used for the Newton Method. - *
- * - * @note if dim - struct TransformR2UInitialGuess - { - static const double KA[GeometryInfo::vertices_per_cell][dim]; - static const double Kb[GeometryInfo::vertices_per_cell]; - }; - - - /* - Octave code: - M=[0 1; 1 1]; - K1 = transpose(M) * inverse (M*transpose(M)); - printf ("{%f, %f},\n", K1' ); - */ - template <> - const double - TransformR2UInitialGuess<1>:: - KA[GeometryInfo<1>::vertices_per_cell][1] = - { - {-1.000000}, - {1.000000} - }; - - template <> - const double - TransformR2UInitialGuess<1>:: - Kb[GeometryInfo<1>::vertices_per_cell] = {1.000000, 0.000000}; - - - /* - Octave code: - M=[0 1 0 1;0 0 1 1;1 1 1 1]; - K2 = transpose(M) * inverse (M*transpose(M)); - printf ("{%f, %f, %f},\n", K2' ); - */ - template <> - const double - TransformR2UInitialGuess<2>:: - KA[GeometryInfo<2>::vertices_per_cell][2] = - { - {-0.500000, -0.500000}, - { 0.500000, -0.500000}, - {-0.500000, 0.500000}, - { 0.500000, 0.500000} - }; - - /* - Octave code: - M=[0 1 0 1 0 1 0 1;0 0 1 1 0 0 1 1; 0 0 0 0 1 1 1 1; 1 1 1 1 1 1 1 1]; - K3 = transpose(M) * inverse (M*transpose(M)) - printf ("{%f, %f, %f, %f},\n", K3' ); - */ - template <> - const double - TransformR2UInitialGuess<2>:: - Kb[GeometryInfo<2>::vertices_per_cell] = - {0.750000,0.250000,0.250000,-0.250000 }; - - - template <> - const double - TransformR2UInitialGuess<3>:: - KA[GeometryInfo<3>::vertices_per_cell][3] = - { - {-0.250000, -0.250000, -0.250000}, - { 0.250000, -0.250000, -0.250000}, - {-0.250000, 0.250000, -0.250000}, - { 0.250000, 0.250000, -0.250000}, - {-0.250000, -0.250000, 0.250000}, - { 0.250000, -0.250000, 0.250000}, - {-0.250000, 0.250000, 0.250000}, - { 0.250000, 0.250000, 0.250000} - - }; - - - template <> - const double - TransformR2UInitialGuess<3>:: - Kb[GeometryInfo<3>::vertices_per_cell] = - {0.500000,0.250000,0.250000,0.000000,0.250000,0.000000,0.000000,-0.250000}; - - template - Point - transform_real_to_unit_cell_initial_guess (const std::vector > &vertex, - const Point &p) - { - Point p_unit; - - dealii::FullMatrix KA(GeometryInfo::vertices_per_cell, dim); - dealii::Vector Kb(GeometryInfo::vertices_per_cell); - - KA.fill( (double *)(TransformR2UInitialGuess::KA) ); - for (unsigned int i=0; i::vertices_per_cell; ++i) - Kb(i) = TransformR2UInitialGuess::Kb[i]; - - FullMatrix Y(spacedim, GeometryInfo::vertices_per_cell); - for (unsigned int v=0; v::vertices_per_cell; v++) - for (unsigned int i=0; i A(spacedim,dim); - Y.mmult(A,KA); // A = Y*KA - dealii::Vector b(spacedim); - Y.vmult(b,Kb); // b = Y*Kb - - for (unsigned int i=0; i dest(dim); - - FullMatrix A_1(dim,spacedim); - if (dim - void - compute_shape_function_values (const unsigned int n_shape_functions, - const std::vector > &unit_points, - typename dealii::MappingManifold<1,spacedim>::InternalData &data) - { - (void)n_shape_functions; - const unsigned int n_points=unit_points.size(); - for (unsigned int k = 0 ; k < n_points ; ++k) - { - double x = unit_points[k](0); - - if (data.shape_values.size()!=0) - { - Assert(data.shape_values.size()==n_shape_functions*n_points, - ExcInternalError()); - data.shape(k,0) = 1.-x; - data.shape(k,1) = x; - } - if (data.shape_derivatives.size()!=0) - { - Assert(data.shape_derivatives.size()==n_shape_functions*n_points, - ExcInternalError()); - data.derivative(k,0)[0] = -1.; - data.derivative(k,1)[0] = 1.; - } - if (data.shape_second_derivatives.size()!=0) - { - // the following may or may not - // work if dim != spacedim - Assert (spacedim == 1, ExcNotImplemented()); - - Assert(data.shape_second_derivatives.size()==n_shape_functions*n_points, - ExcInternalError()); - data.second_derivative(k,0)[0][0] = 0; - data.second_derivative(k,1)[0][0] = 0; - } - if (data.shape_third_derivatives.size()!=0) - { - // if lower order derivative don't work, neither should this - Assert (spacedim == 1, ExcNotImplemented()); - - Assert(data.shape_third_derivatives.size()==n_shape_functions*n_points, - ExcInternalError()); - - Tensor<3,1> zero; - data.third_derivative(k,0) = zero; - data.third_derivative(k,1) = zero; - } - if (data.shape_fourth_derivatives.size()!=0) - { - // if lower order derivative don't work, neither should this - Assert (spacedim == 1, ExcNotImplemented()); - - Assert(data.shape_fourth_derivatives.size()==n_shape_functions*n_points, - ExcInternalError()); - - Tensor<4,1> zero; - data.fourth_derivative(k,0) = zero; - data.fourth_derivative(k,1) = zero; - } - } - } - - - template - void - compute_shape_function_values (const unsigned int n_shape_functions, - const std::vector > &unit_points, - typename dealii::MappingManifold<2,spacedim>::InternalData &data) - { - (void)n_shape_functions; - const unsigned int n_points=unit_points.size(); - for (unsigned int k = 0 ; k < n_points ; ++k) - { - double x = unit_points[k](0); - double y = unit_points[k](1); - - if (data.shape_values.size()!=0) - { - Assert(data.shape_values.size()==n_shape_functions*n_points, - ExcInternalError()); - data.shape(k,0) = (1.-x)*(1.-y); - data.shape(k,1) = x*(1.-y); - data.shape(k,2) = (1.-x)*y; - data.shape(k,3) = x*y; - } - if (data.shape_derivatives.size()!=0) - { - Assert(data.shape_derivatives.size()==n_shape_functions*n_points, - ExcInternalError()); - data.derivative(k,0)[0] = (y-1.); - data.derivative(k,1)[0] = (1.-y); - data.derivative(k,2)[0] = -y; - data.derivative(k,3)[0] = y; - data.derivative(k,0)[1] = (x-1.); - data.derivative(k,1)[1] = -x; - data.derivative(k,2)[1] = (1.-x); - data.derivative(k,3)[1] = x; - } - if (data.shape_second_derivatives.size()!=0) - { - Assert(data.shape_second_derivatives.size()==n_shape_functions*n_points, - ExcInternalError()); - data.second_derivative(k,0)[0][0] = 0; - data.second_derivative(k,1)[0][0] = 0; - data.second_derivative(k,2)[0][0] = 0; - data.second_derivative(k,3)[0][0] = 0; - data.second_derivative(k,0)[0][1] = 1.; - data.second_derivative(k,1)[0][1] = -1.; - data.second_derivative(k,2)[0][1] = -1.; - data.second_derivative(k,3)[0][1] = 1.; - data.second_derivative(k,0)[1][0] = 1.; - data.second_derivative(k,1)[1][0] = -1.; - data.second_derivative(k,2)[1][0] = -1.; - data.second_derivative(k,3)[1][0] = 1.; - data.second_derivative(k,0)[1][1] = 0; - data.second_derivative(k,1)[1][1] = 0; - data.second_derivative(k,2)[1][1] = 0; - data.second_derivative(k,3)[1][1] = 0; - } - if (data.shape_third_derivatives.size()!=0) - { - Assert(data.shape_third_derivatives.size()==n_shape_functions*n_points, - ExcInternalError()); - - Tensor<3,2> zero; - for (unsigned int i=0; i<4; ++i) - data.third_derivative(k,i) = zero; - } - if (data.shape_fourth_derivatives.size()!=0) - { - Assert(data.shape_fourth_derivatives.size()==n_shape_functions*n_points, - ExcInternalError()); - Tensor<4,2> zero; - for (unsigned int i=0; i<4; ++i) - data.fourth_derivative(k,i) = zero; - } - } - } - - - - template - void - compute_shape_function_values (const unsigned int n_shape_functions, - const std::vector > &unit_points, - typename dealii::MappingManifold<3,spacedim>::InternalData &data) - { - (void)n_shape_functions; - const unsigned int n_points=unit_points.size(); - for (unsigned int k = 0 ; k < n_points ; ++k) - { - double x = unit_points[k](0); - double y = unit_points[k](1); - double z = unit_points[k](2); - - if (data.shape_values.size()!=0) - { - Assert(data.shape_values.size()==n_shape_functions*n_points, - ExcInternalError()); - data.shape(k,0) = (1.-x)*(1.-y)*(1.-z); - data.shape(k,1) = x*(1.-y)*(1.-z); - data.shape(k,2) = (1.-x)*y*(1.-z); - data.shape(k,3) = x*y*(1.-z); - data.shape(k,4) = (1.-x)*(1.-y)*z; - data.shape(k,5) = x*(1.-y)*z; - data.shape(k,6) = (1.-x)*y*z; - data.shape(k,7) = x*y*z; - } - if (data.shape_derivatives.size()!=0) - { - Assert(data.shape_derivatives.size()==n_shape_functions*n_points, - ExcInternalError()); - data.derivative(k,0)[0] = (y-1.)*(1.-z); - data.derivative(k,1)[0] = (1.-y)*(1.-z); - data.derivative(k,2)[0] = -y*(1.-z); - data.derivative(k,3)[0] = y*(1.-z); - data.derivative(k,4)[0] = (y-1.)*z; - data.derivative(k,5)[0] = (1.-y)*z; - data.derivative(k,6)[0] = -y*z; - data.derivative(k,7)[0] = y*z; - data.derivative(k,0)[1] = (x-1.)*(1.-z); - data.derivative(k,1)[1] = -x*(1.-z); - data.derivative(k,2)[1] = (1.-x)*(1.-z); - data.derivative(k,3)[1] = x*(1.-z); - data.derivative(k,4)[1] = (x-1.)*z; - data.derivative(k,5)[1] = -x*z; - data.derivative(k,6)[1] = (1.-x)*z; - data.derivative(k,7)[1] = x*z; - data.derivative(k,0)[2] = (x-1)*(1.-y); - data.derivative(k,1)[2] = x*(y-1.); - data.derivative(k,2)[2] = (x-1.)*y; - data.derivative(k,3)[2] = -x*y; - data.derivative(k,4)[2] = (1.-x)*(1.-y); - data.derivative(k,5)[2] = x*(1.-y); - data.derivative(k,6)[2] = (1.-x)*y; - data.derivative(k,7)[2] = x*y; - } - if (data.shape_second_derivatives.size()!=0) - { - // the following may or may not - // work if dim != spacedim - Assert (spacedim == 3, ExcNotImplemented()); - - Assert(data.shape_second_derivatives.size()==n_shape_functions*n_points, - ExcInternalError()); - data.second_derivative(k,0)[0][0] = 0; - data.second_derivative(k,1)[0][0] = 0; - data.second_derivative(k,2)[0][0] = 0; - data.second_derivative(k,3)[0][0] = 0; - data.second_derivative(k,4)[0][0] = 0; - data.second_derivative(k,5)[0][0] = 0; - data.second_derivative(k,6)[0][0] = 0; - data.second_derivative(k,7)[0][0] = 0; - data.second_derivative(k,0)[1][1] = 0; - data.second_derivative(k,1)[1][1] = 0; - data.second_derivative(k,2)[1][1] = 0; - data.second_derivative(k,3)[1][1] = 0; - data.second_derivative(k,4)[1][1] = 0; - data.second_derivative(k,5)[1][1] = 0; - data.second_derivative(k,6)[1][1] = 0; - data.second_derivative(k,7)[1][1] = 0; - data.second_derivative(k,0)[2][2] = 0; - data.second_derivative(k,1)[2][2] = 0; - data.second_derivative(k,2)[2][2] = 0; - data.second_derivative(k,3)[2][2] = 0; - data.second_derivative(k,4)[2][2] = 0; - data.second_derivative(k,5)[2][2] = 0; - data.second_derivative(k,6)[2][2] = 0; - data.second_derivative(k,7)[2][2] = 0; - - data.second_derivative(k,0)[0][1] = (1.-z); - data.second_derivative(k,1)[0][1] = -(1.-z); - data.second_derivative(k,2)[0][1] = -(1.-z); - data.second_derivative(k,3)[0][1] = (1.-z); - data.second_derivative(k,4)[0][1] = z; - data.second_derivative(k,5)[0][1] = -z; - data.second_derivative(k,6)[0][1] = -z; - data.second_derivative(k,7)[0][1] = z; - data.second_derivative(k,0)[1][0] = (1.-z); - data.second_derivative(k,1)[1][0] = -(1.-z); - data.second_derivative(k,2)[1][0] = -(1.-z); - data.second_derivative(k,3)[1][0] = (1.-z); - data.second_derivative(k,4)[1][0] = z; - data.second_derivative(k,5)[1][0] = -z; - data.second_derivative(k,6)[1][0] = -z; - data.second_derivative(k,7)[1][0] = z; - - data.second_derivative(k,0)[0][2] = (1.-y); - data.second_derivative(k,1)[0][2] = -(1.-y); - data.second_derivative(k,2)[0][2] = y; - data.second_derivative(k,3)[0][2] = -y; - data.second_derivative(k,4)[0][2] = -(1.-y); - data.second_derivative(k,5)[0][2] = (1.-y); - data.second_derivative(k,6)[0][2] = -y; - data.second_derivative(k,7)[0][2] = y; - data.second_derivative(k,0)[2][0] = (1.-y); - data.second_derivative(k,1)[2][0] = -(1.-y); - data.second_derivative(k,2)[2][0] = y; - data.second_derivative(k,3)[2][0] = -y; - data.second_derivative(k,4)[2][0] = -(1.-y); - data.second_derivative(k,5)[2][0] = (1.-y); - data.second_derivative(k,6)[2][0] = -y; - data.second_derivative(k,7)[2][0] = y; - - data.second_derivative(k,0)[1][2] = (1.-x); - data.second_derivative(k,1)[1][2] = x; - data.second_derivative(k,2)[1][2] = -(1.-x); - data.second_derivative(k,3)[1][2] = -x; - data.second_derivative(k,4)[1][2] = -(1.-x); - data.second_derivative(k,5)[1][2] = -x; - data.second_derivative(k,6)[1][2] = (1.-x); - data.second_derivative(k,7)[1][2] = x; - data.second_derivative(k,0)[2][1] = (1.-x); - data.second_derivative(k,1)[2][1] = x; - data.second_derivative(k,2)[2][1] = -(1.-x); - data.second_derivative(k,3)[2][1] = -x; - data.second_derivative(k,4)[2][1] = -(1.-x); - data.second_derivative(k,5)[2][1] = -x; - data.second_derivative(k,6)[2][1] = (1.-x); - data.second_derivative(k,7)[2][1] = x; - } - if (data.shape_third_derivatives.size()!=0) - { - // if lower order derivative don't work, neither should this - Assert (spacedim == 3, ExcNotImplemented()); - - Assert(data.shape_third_derivatives.size()==n_shape_functions*n_points, - ExcInternalError()); - - for (unsigned int i=0; i<3; ++i) - for (unsigned int j=0; j<3; ++j) - for (unsigned int l=0; l<3; ++l) - if ((i==j)||(j==l)||(l==i)) - { - for (unsigned int m=0; m<8; ++m) - data.third_derivative(k,m)[i][j][l] = 0; - } - else - { - data.third_derivative(k,0)[i][j][l] = -1.; - data.third_derivative(k,1)[i][j][l] = 1.; - data.third_derivative(k,2)[i][j][l] = 1.; - data.third_derivative(k,3)[i][j][l] = -1.; - data.third_derivative(k,4)[i][j][l] = 1.; - data.third_derivative(k,5)[i][j][l] = -1.; - data.third_derivative(k,6)[i][j][l] = -1.; - data.third_derivative(k,7)[i][j][l] = 1.; - } - - } - if (data.shape_fourth_derivatives.size()!=0) - { - // if lower order derivative don't work, neither should this - Assert (spacedim == 3, ExcNotImplemented()); - - Assert(data.shape_fourth_derivatives.size()==n_shape_functions*n_points, - ExcInternalError()); - Tensor<4,3> zero; - for (unsigned int i=0; i<8; ++i) - data.fourth_derivative(k,i) = zero; - } - } - } - } - } -} - - - - - template MappingManifold::InternalData::InternalData () - : - polynomial_degree (1), - n_shape_functions (Utilities::fixed_power(polynomial_degree+1)) {} @@ -674,18 +50,18 @@ template std::size_t MappingManifold::InternalData::memory_consumption () const { - return (Mapping::InternalDataBase::memory_consumption() + - MemoryConsumption::memory_consumption (shape_values) + - MemoryConsumption::memory_consumption (shape_derivatives) + - MemoryConsumption::memory_consumption (covariant) + - MemoryConsumption::memory_consumption (contravariant) + - MemoryConsumption::memory_consumption (unit_tangentials) + - MemoryConsumption::memory_consumption (aux) + - MemoryConsumption::memory_consumption (mapping_support_points) + - MemoryConsumption::memory_consumption (cell_of_current_support_points) + - MemoryConsumption::memory_consumption (volume_elements) + - MemoryConsumption::memory_consumption (polynomial_degree) + - MemoryConsumption::memory_consumption (n_shape_functions)); + return (Mapping::InternalDataBase::memory_consumption() ); + // MemoryConsumption::memory_consumption (shape_values) + + // MemoryConsumption::memory_consumption (shape_derivatives) + + // MemoryConsumption::memory_consumption (covariant) + + // MemoryConsumption::memory_consumption (contravariant) + + // MemoryConsumption::memory_consumption (unit_tangentials) + + // MemoryConsumption::memory_consumption (aux) + + // MemoryConsumption::memory_consumption (mapping_support_points) + + // MemoryConsumption::memory_consumption (cell_of_current_support_points) + + // MemoryConsumption::memory_consumption (volume_elements) + + // MemoryConsumption::memory_consumption (polynomial_degree) + + // MemoryConsumption::memory_consumption (n_shape_functions)); } @@ -702,49 +78,52 @@ initialize (const UpdateFlags update_flags, const unsigned int n_q_points = q.size(); + // Update the weights used in the Manifold Quadrature formulas + compute_manifold_quadrature_weights(q); + // see if we need the (transformation) shape function values // and/or gradients and resize the necessary arrays if (this->update_each & update_quadrature_points) - shape_values.resize(n_shape_functions * n_q_points); - - if (this->update_each & (update_covariant_transformation - | update_contravariant_transformation - | update_JxW_values - | update_boundary_forms - | update_normal_vectors - | update_jacobians - | update_jacobian_grads - | update_inverse_jacobians - | update_jacobian_pushed_forward_grads - | update_jacobian_2nd_derivatives - | update_jacobian_pushed_forward_2nd_derivatives - | update_jacobian_3rd_derivatives - | update_jacobian_pushed_forward_3rd_derivatives)) - shape_derivatives.resize(n_shape_functions * n_q_points); - - if (this->update_each & update_covariant_transformation) - covariant.resize(n_original_q_points); - - if (this->update_each & update_contravariant_transformation) - contravariant.resize(n_original_q_points); - - if (this->update_each & update_volume_elements) - volume_elements.resize(n_original_q_points); - - if (this->update_each & - (update_jacobian_grads | update_jacobian_pushed_forward_grads) ) - shape_second_derivatives.resize(n_shape_functions * n_q_points); - - if (this->update_each & - (update_jacobian_2nd_derivatives | update_jacobian_pushed_forward_2nd_derivatives) ) - shape_third_derivatives.resize(n_shape_functions * n_q_points); - - if (this->update_each & - (update_jacobian_3rd_derivatives | update_jacobian_pushed_forward_3rd_derivatives) ) - shape_fourth_derivatives.resize(n_shape_functions * n_q_points); - - // now also fill the various fields with their correct values - compute_shape_function_values (q.get_points()); + cell_manifold_quadratures.resize(q.size()); + + // if (this->update_each & (update_covariant_transformation + // | update_contravariant_transformation + // | update_JxW_values + // | update_boundary_forms + // | update_normal_vectors + // | update_jacobians + // | update_jacobian_grads + // | update_inverse_jacobians + // | update_jacobian_pushed_forward_grads + // | update_jacobian_2nd_derivatives + // | update_jacobian_pushed_forward_2nd_derivatives + // | update_jacobian_3rd_derivatives + // | update_jacobian_pushed_forward_3rd_derivatives)) + // shape_derivatives.resize(n_shape_functions * n_q_points); + + // if (this->update_each & update_covariant_transformation) + // covariant.resize(n_original_q_points); + + // if (this->update_each & update_contravariant_transformation) + // contravariant.resize(n_original_q_points); + + // if (this->update_each & update_volume_elements) + // volume_elements.resize(n_original_q_points); + + // if (this->update_each & + // (update_jacobian_grads | update_jacobian_pushed_forward_grads) ) + // shape_second_derivatives.resize(n_shape_functions * n_q_points); + + // if (this->update_each & + // (update_jacobian_2nd_derivatives | update_jacobian_pushed_forward_2nd_derivatives) ) + // shape_third_derivatives.resize(n_shape_functions * n_q_points); + + // if (this->update_each & + // (update_jacobian_3rd_derivatives | update_jacobian_pushed_forward_3rd_derivatives) ) + // shape_fourth_derivatives.resize(n_shape_functions * n_q_points); + + // // now also fill the various fields with their correct values + // compute_shape_function_values (q.get_points()); } @@ -756,413 +135,72 @@ initialize_face (const UpdateFlags update_flags, const Quadrature &q, const unsigned int n_original_q_points) { - initialize (update_flags, q, n_original_q_points); - - if (dim > 1) - { - if (this->update_each & update_boundary_forms) - { - aux.resize (dim-1, std::vector > (n_original_q_points)); - - // Compute tangentials to the - // unit cell. - const unsigned int nfaces = GeometryInfo::faces_per_cell; - unit_tangentials.resize (nfaces*(dim-1), - std::vector > (n_original_q_points)); - if (dim==2) - { - // ensure a counterclockwise - // orientation of tangentials - static const int tangential_orientation[4]= {-1,1,1,-1}; - for (unsigned int i=0; i tang; - tang[1-i/2]=tangential_orientation[i]; - std::fill (unit_tangentials[i].begin(), - unit_tangentials[i].end(), tang); - } - } - else if (dim==3) - { - for (unsigned int i=0; i tang1, tang2; - - const unsigned int nd= - GeometryInfo::unit_normal_direction[i]; - - // first tangential - // vector in direction - // of the (nd+1)%3 axis - // and inverted in case - // of unit inward normal - tang1[(nd+1)%dim]=GeometryInfo::unit_normal_orientation[i]; - // second tangential - // vector in direction - // of the (nd+2)%3 axis - tang2[(nd+2)%dim]=1.; - - // same unit tangents - // for all quadrature - // points on this face - std::fill (unit_tangentials[i].begin(), - unit_tangentials[i].end(), tang1); - std::fill (unit_tangentials[nfaces+i].begin(), - unit_tangentials[nfaces+i].end(), tang2); - } - } - } - } -} - - - -namespace -{ - template - std::vector - get_dpo_vector () - { - unsigned int degree = 1; - std::vector dpo(dim+1, 1U); - for (unsigned int i=1; i 1) + // { + // if (this->update_each & update_boundary_forms) + // { + // aux.resize (dim-1, std::vector > (n_original_q_points)); + + // // Compute tangentials to the + // // unit cell. + // const unsigned int nfaces = GeometryInfo::faces_per_cell; + // unit_tangentials.resize (nfaces*(dim-1), + // std::vector > (n_original_q_points)); + // if (dim==2) + // { + // // ensure a counterclockwise + // // orientation of tangentials + // static const int tangential_orientation[4]= {-1,1,1,-1}; + // for (unsigned int i=0; i tang; + // tang[1-i/2]=tangential_orientation[i]; + // std::fill (unit_tangentials[i].begin(), + // unit_tangentials[i].end(), tang); + // } + // } + // else if (dim==3) + // { + // for (unsigned int i=0; i tang1, tang2; + + // const unsigned int nd= + // GeometryInfo::unit_normal_direction[i]; + + // // first tangential + // // vector in direction + // // of the (nd+1)%3 axis + // // and inverted in case + // // of unit inward normal + // tang1[(nd+1)%dim]=GeometryInfo::unit_normal_orientation[i]; + // // second tangential + // // vector in direction + // // of the (nd+2)%3 axis + // tang2[(nd+2)%dim]=1.; + + // // same unit tangents + // // for all quadrature + // // points on this face + // std::fill (unit_tangentials[i].begin(), + // unit_tangentials[i].end(), tang1); + // std::fill (unit_tangentials[nfaces+i].begin(), + // unit_tangentials[nfaces+i].end(), tang2); + // } + // } + // } + // } } - -template -void -MappingManifold::InternalData:: -compute_shape_function_values (const std::vector > &unit_points) -{ - // if the polynomial degree is one, then we can simplify code a bit - // by using hard-coded shape functions. - if ((polynomial_degree == 1) - && - (dim == spacedim)) - internal::MappingQ1::compute_shape_function_values (n_shape_functions, - unit_points, *this); - else - // otherwise ask an object that describes the polynomial space - { - const unsigned int n_points=unit_points.size(); - - // Construct the tensor product polynomials used as shape functions for the - // Qp mapping of cells at the boundary. - const QGaussLobatto<1> line_support_points (polynomial_degree + 1); - const TensorProductPolynomials - tensor_pols (Polynomials::generate_complete_Lagrange_basis(line_support_points.get_points())); - Assert (n_shape_functions==tensor_pols.n(), - ExcInternalError()); - - // then also construct the mapping from lexicographic to the Qp shape function numbering - const std::vector - renumber (FETools:: - lexicographic_to_hierarchic_numbering ( - FiniteElementData (get_dpo_vector(), 1, - polynomial_degree))); - - std::vector values; - std::vector > grads; - if (shape_values.size()!=0) - { - Assert(shape_values.size()==n_shape_functions*n_points, - ExcInternalError()); - values.resize(n_shape_functions); - } - if (shape_derivatives.size()!=0) - { - Assert(shape_derivatives.size()==n_shape_functions*n_points, - ExcInternalError()); - grads.resize(n_shape_functions); - } - - std::vector > grad2; - if (shape_second_derivatives.size()!=0) - { - Assert(shape_second_derivatives.size()==n_shape_functions*n_points, - ExcInternalError()); - grad2.resize(n_shape_functions); - } - - std::vector > grad3; - if (shape_third_derivatives.size()!=0) - { - Assert(shape_third_derivatives.size()==n_shape_functions*n_points, - ExcInternalError()); - grad3.resize(n_shape_functions); - } - - std::vector > grad4; - if (shape_fourth_derivatives.size()!=0) - { - Assert(shape_fourth_derivatives.size()==n_shape_functions*n_points, - ExcInternalError()); - grad4.resize(n_shape_functions); - } - - - if (shape_values.size()!=0 || - shape_derivatives.size()!=0 || - shape_second_derivatives.size()!=0 || - shape_third_derivatives.size()!=0 || - shape_fourth_derivatives.size()!=0 ) - for (unsigned int point=0; pointsupport_point_weights_on_quad(hex) arrays. - * - * Called by the compute_support_point_weights_on_quad(hex) functions if the - * data is not yet hardcoded. - * - * For the definition of the support_point_weights_on_quad(hex) please - * refer to equation (8) of the `mapping' report. - */ - template - Table<2,double> - compute_laplace_vector(const unsigned int polynomial_degree) - { - Table<2,double> lvs; - - Assert(lvs.n_rows()==0, ExcInternalError()); - Assert(dim==2 || dim==3, ExcNotImplemented()); - - // for degree==1, we shouldn't have to compute any support points, since all - // of them are on the vertices - Assert(polynomial_degree>1, ExcInternalError()); - - const unsigned int n_inner = Utilities::fixed_power(polynomial_degree-1); - const unsigned int n_outer = (dim==1) ? 2 : - ((dim==2) ? - 4+4*(polynomial_degree-1) : - 8+12*(polynomial_degree-1)+6*(polynomial_degree-1)*(polynomial_degree-1)); - - - // compute the shape gradients at the quadrature points on the unit cell - const QGauss quadrature(polynomial_degree+1); - const unsigned int n_q_points=quadrature.size(); - - typename MappingManifold::InternalData quadrature_data(polynomial_degree); - quadrature_data.shape_derivatives.resize(quadrature_data.n_shape_functions * - n_q_points); - quadrature_data.compute_shape_function_values(quadrature.get_points()); - - // Compute the stiffness matrix of the inner dofs - FullMatrix S(n_inner); - for (unsigned int point=0; point T(n_inner, n_outer); - for (unsigned int point=0; point S_1(n_inner); - S_1.invert(S); - - FullMatrix S_1_T(n_inner, n_outer); - - // S:=S_1*T - S_1.mmult(S_1_T,T); - - // Resize and initialize the lvs - lvs.reinit (n_inner, n_outer); - for (unsigned int i=0; iMappingQ for dim= 2 and 3. - * - * For degree<4 this function sets the @p support_point_weights_on_quad to - * the hardcoded data. For degree>=4 and MappingQ<2> this vector is - * computed. - * - * For the definition of the @p support_point_weights_on_quad please refer to - * equation (8) of the `mapping' report. - */ - template - Table<2,double> - compute_support_point_weights_on_quad(const unsigned int polynomial_degree) - { - Table<2,double> loqvs; - - // in 1d, there are no quads, so return an empty object - if (dim == 1) - return loqvs; - - // we are asked to compute weights for interior support points, but - // there are no interior points if degree==1 - if (polynomial_degree == 1) - return loqvs; - - const unsigned int n_inner_2d=(polynomial_degree-1)*(polynomial_degree-1); - const unsigned int n_outer_2d=4+4*(polynomial_degree-1); - - // first check whether we have precomputed the values for some polynomial - // degree; the sizes of arrays is n_inner_2d*n_outer_2d - if (polynomial_degree == 2) - { - // (checked these values against the output of compute_laplace_vector - // again, and found they're indeed right -- just in case someone wonders - // where they come from -- WB) - static const double loqv2[1*8] - = {1/16., 1/16., 1/16., 1/16., 3/16., 3/16., 3/16., 3/16.}; - Assert (sizeof(loqv2)/sizeof(loqv2[0]) == - n_inner_2d * n_outer_2d, - ExcInternalError()); - - // copy and return - loqvs.reinit(n_inner_2d, n_outer_2d); - for (unsigned int unit_point=0; unit_point(polynomial_degree); - } - - // the sum of weights of the points at the outer rim should be one. check - // this - for (unsigned int unit_point=0; unit_pointMappingQ<3>. - * - * For degree==2 this function sets the @p support_point_weights_on_hex to - * the hardcoded data. For degree>2 this vector is computed. - * - * For the definition of the @p support_point_weights_on_hex please refer to - * equation (8) of the `mapping' report. - */ - template - Table<2,double> - compute_support_point_weights_on_hex(const unsigned int polynomial_degree) - { - Table<2,double> lohvs; - - // in 1d and 2d, there are no hexes, so return an empty object - if (dim < 3) - return lohvs; - - // we are asked to compute weights for interior support points, but - // there are no interior points if degree==1 - if (polynomial_degree == 1) - return lohvs; - - const unsigned int n_inner = Utilities::fixed_power(polynomial_degree-1); - const unsigned int n_outer = 8+12*(polynomial_degree-1)+6*(polynomial_degree-1)*(polynomial_degree-1); - - // first check whether we have precomputed the values for some polynomial - // degree; the sizes of arrays is n_inner_2d*n_outer_2d - if (polynomial_degree == 2) - { - static const double lohv2[26] - = {1/128., 1/128., 1/128., 1/128., 1/128., 1/128., 1/128., 1/128., - 7/192., 7/192., 7/192., 7/192., 7/192., 7/192., 7/192., 7/192., - 7/192., 7/192., 7/192., 7/192., - 1/12., 1/12., 1/12., 1/12., 1/12., 1/12. - }; - - // copy and return - lohvs.reinit(n_inner, n_outer); - for (unsigned int unit_point=0; unit_point(polynomial_degree); - } - - // the sum of weights of the points at the outer rim should be one. check - // this - for (unsigned int unit_point=0; unit_point MappingManifold::MappingManifold () : - polynomial_degree(1), - line_support_points(this->polynomial_degree+1), - fe_q(dim == 3 ? new FE_Q(this->polynomial_degree) : 0), - support_point_weights_on_quad (compute_support_point_weights_on_quad(this->polynomial_degree)), - support_point_weights_on_hex (compute_support_point_weights_on_hex(this->polynomial_degree)) + fe_q(1) + // support_point_weights_on_quad (compute_support_point_weights_on_quad(this->polynomial_degree)), + // support_point_weights_on_hex (compute_support_point_weights_on_hex(this->polynomial_degree)), { } @@ -1171,11 +209,7 @@ MappingManifold::MappingManifold () template MappingManifold::MappingManifold (const MappingManifold &mapping) : - polynomial_degree(mapping.polynomial_degree), - line_support_points(mapping.line_support_points), - fe_q(dim == 3 ? new FE_Q(*mapping.fe_q) : 0), - support_point_weights_on_quad (mapping.support_point_weights_on_quad), - support_point_weights_on_hex (mapping.support_point_weights_on_hex) + fe_q(1) {} @@ -1196,28 +230,14 @@ MappingManifold:: transform_unit_to_real_cell (const typename Triangulation::cell_iterator &cell, const Point &p) const { - // set up the polynomial space - const QGaussLobatto<1> line_support_points (polynomial_degree + 1); - const TensorProductPolynomials - tensor_pols (Polynomials::generate_complete_Lagrange_basis(line_support_points.get_points())); - Assert (tensor_pols.n() == Utilities::fixed_power(polynomial_degree+1), - ExcInternalError()); - - // then also construct the mapping from lexicographic to the Qp shape function numbering - const std::vector - renumber (FETools:: - lexicographic_to_hierarchic_numbering ( - FiniteElementData (get_dpo_vector(polynomial_degree), 1, - polynomial_degree))); - - const std::vector > support_points - = this->compute_mapping_support_points(cell); - - Point mapped_point; - for (unsigned int i=0; i > vertices; + std::vector > weights; + for (unsigned int v=0; v::vertices_per_cell; ++v) + { + vertices.push_back(cell->vertex(v)); + weights.push_back(fe_q.shape_value(v,p)); + } + return cell->get_manifold().get_new_point(Quadrature(vertices, weights)); } @@ -1232,629 +252,6 @@ template <> class MappingManifold<3,4> {}; -namespace -{ - /** - * Using the relative weights of the shape functions evaluated at - * one point on the reference cell (and stored in data.shape_values - * and accessed via data.shape(0,i)) and the locations of mapping - * support points (stored in data.mapping_support_points), compute - * the mapped location of that point in real space. - */ - template - Point - compute_mapped_location_of_point (const typename MappingManifold::InternalData &data) - { - AssertDimension (data.shape_values.size(), - data.mapping_support_points.size()); - - // use now the InternalData to compute the point in real space. - Point p_real; - for (unsigned int i=0; i - Point - do_transform_real_to_unit_cell_internal - (const typename Triangulation::cell_iterator &cell, - const Point &p, - const Point &initial_p_unit, - typename MappingManifold::InternalData &mdata) - { - const unsigned int spacedim = dim; - - const unsigned int n_shapes=mdata.shape_values.size(); - (void)n_shapes; - Assert(n_shapes!=0, ExcInternalError()); - AssertDimension (mdata.shape_derivatives.size(), n_shapes); - - std::vector > &points=mdata.mapping_support_points; - AssertDimension (points.size(), n_shapes); - - - // Newton iteration to solve - // f(x)=p(x)-p=0 - // where we are looking for 'x' and p(x) is the forward transformation - // from unit to real cell. We solve this using a Newton iteration - // x_{n+1}=x_n-[f'(x)]^{-1}f(x) - // The start value is set to be the linear approximation to the cell - - // The shape values and derivatives of the mapping at this point are - // previously computed. - - Point p_unit = initial_p_unit; - - mdata.compute_shape_function_values(std::vector > (1, p_unit)); - - Point p_real = compute_mapped_location_of_point(mdata); - Tensor<1,spacedim> f = p_real-p; - - // early out if we already have our point - if (f.norm_square() < 1e-24 * cell->diameter() * cell->diameter()) - return p_unit; - - // we need to compare the position of the computed p(x) against the given - // point 'p'. We will terminate the iteration and return 'x' if they are - // less than eps apart. The question is how to choose eps -- or, put maybe - // more generally: in which norm we want these 'p' and 'p(x)' to be eps - // apart. - // - // the question is difficult since we may have to deal with very elongated - // cells where we may achieve 1e-12*h for the distance of these two points - // in the 'long' direction, but achieving this tolerance in the 'short' - // direction of the cell may not be possible - // - // what we do instead is then to terminate iterations if - // \| p(x) - p \|_A < eps - // where the A-norm is somehow induced by the transformation of the cell. - // in particular, we want to measure distances relative to the sizes of - // the cell in its principal directions. - // - // to define what exactly A should be, note that to first order we have - // the following (assuming that x* is the solution of the problem, i.e., - // p(x*)=p): - // p(x) - p = p(x) - p(x*) - // = -grad p(x) * (x*-x) + higher order terms - // This suggest to measure with a norm that corresponds to - // A = {[grad p(x]^T [grad p(x)]}^{-1} - // because then - // \| p(x) - p \|_A \approx \| x - x* \| - // Consequently, we will try to enforce that - // \| p(x) - p \|_A = \| f \| <= eps - // - // Note that using this norm is a bit dangerous since the norm changes - // in every iteration (A isn't fixed by depends on xk). However, if the - // cell is not too deformed (it may be stretched, but not twisted) then - // the mapping is almost linear and A is indeed constant or nearly so. - const double eps = 1.e-11; - const unsigned int newton_iteration_limit = 20; - - unsigned int newton_iteration = 0; - double last_f_weighted_norm; - do - { -#ifdef DEBUG_TRANSFORM_REAL_TO_UNIT_CELL - std::cout << "Newton iteration " << newton_iteration << std::endl; -#endif - - // f'(x) - Tensor<2,spacedim> df; - for (unsigned int k=0; k &grad_transform=mdata.derivative(0,k); - const Point &point=points[k]; - - for (unsigned int i=0; i df_inverse = invert(df); - const Tensor<1,spacedim> delta = df_inverse * static_cast&>(f); - -#ifdef DEBUG_TRANSFORM_REAL_TO_UNIT_CELL - std::cout << " delta=" << delta << std::endl; -#endif - - // do a line search - double step_length = 1; - do - { - // update of p_unit. The spacedim-th component of transformed point - // is simply ignored in codimension one case. When this component is - // not zero, then we are projecting the point to the surface or - // curve identified by the cell. - Point p_unit_trial = p_unit; - for (unsigned int i=0; i > (1, p_unit_trial)); - - // f(x) - Point p_real_trial = compute_mapped_location_of_point(mdata); - const Tensor<1,spacedim> f_trial = p_real_trial-p; - -#ifdef DEBUG_TRANSFORM_REAL_TO_UNIT_CELL - std::cout << " step_length=" << step_length << std::endl - << " ||f || =" << f.norm() << std::endl - << " ||f*|| =" << f_trial.norm() << std::endl - << " ||f*||_A =" << (df_inverse * f_trial).norm() << std::endl; -#endif - - // see if we are making progress with the current step length - // and if not, reduce it by a factor of two and try again - // - // strictly speaking, we should probably use the same norm as we use - // for the outer algorithm. in practice, line search is just a - // crutch to find a "reasonable" step length, and so using the l2 - // norm is probably just fine - if (f_trial.norm() < f.norm()) - { - p_real = p_real_trial; - p_unit = p_unit_trial; - f = f_trial; - break; - } - else if (step_length > 0.05) - step_length /= 2; - else - AssertThrow (false, - (typename Mapping::ExcTransformationFailed())); - } - while (true); - - ++newton_iteration; - if (newton_iteration > newton_iteration_limit) - AssertThrow (false, - (typename Mapping::ExcTransformationFailed())); - last_f_weighted_norm = (df_inverse * f).norm(); - } - while (last_f_weighted_norm > eps); - - return p_unit; - } - - - - /** - * Implementation of transform_real_to_unit_cell for dim==spacedim-1 - */ - template - Point - do_transform_real_to_unit_cell_internal_codim1 - (const typename Triangulation::cell_iterator &cell, - const Point &p, - const Point &initial_p_unit, - typename MappingManifold::InternalData &mdata) - { - const unsigned int spacedim = dim+1; - - const unsigned int n_shapes=mdata.shape_values.size(); - (void)n_shapes; - Assert(n_shapes!=0, ExcInternalError()); - Assert(mdata.shape_derivatives.size()==n_shapes, ExcInternalError()); - Assert(mdata.shape_second_derivatives.size()==n_shapes, ExcInternalError()); - - std::vector > &points=mdata.mapping_support_points; - Assert(points.size()==n_shapes, ExcInternalError()); - - Point p_minus_F; - - Tensor<1,spacedim> DF[dim]; - Tensor<1,spacedim> D2F[dim][dim]; - - Point p_unit = initial_p_unit; - Point f; - Tensor<2,dim> df; - - // Evaluate first and second derivatives - mdata.compute_shape_function_values(std::vector > (1, p_unit)); - - for (unsigned int k=0; k &grad_phi_k = mdata.derivative(0,k); - const Tensor<2,dim> &hessian_k = mdata.second_derivative(0,k); - const Point &point_k = points[k]; - - for (unsigned int j=0; j(mdata); - - - for (unsigned int j=0; jdiameter(); - const unsigned int loop_limit = 10; - - unsigned int loop=0; - - while (f.norm()>eps && loop++ d = invert(df) * static_cast&>(f); - p_unit -= d; - - for (unsigned int j=0; j > (1, p_unit)); - - for (unsigned int k=0; k &grad_phi_k = mdata.derivative(0,k); - const Tensor<2,dim> &hessian_k = mdata.second_derivative(0,k); - const Point &point_k = points[k]; - - for (unsigned int j=0; j(mdata); - - for (unsigned int j=0; j::ExcTransformationFailed())); - - return p_unit; - } - - -} - - - -// visual studio freaks out when trying to determine if -// do_transform_real_to_unit_cell_internal with dim=3 and spacedim=4 is a good -// candidate. So instead of letting the compiler pick the correct overload, we -// use template specialization to make sure we pick up the right function to -// call: - -template -Point -MappingManifold:: -transform_real_to_unit_cell_internal -(const typename Triangulation::cell_iterator &, - const Point &, - const Point &) const -{ - // default implementation (should never be called) - Assert(false, ExcInternalError()); - return Point(); -} - -template<> -Point<1> -MappingManifold<1,1>:: -transform_real_to_unit_cell_internal -(const Triangulation<1,1>::cell_iterator &cell, - const Point<1> &p, - const Point<1> &initial_p_unit) const -{ - const int dim = 1; - const int spacedim = 1; - - const Quadrature point_quadrature(initial_p_unit); - - UpdateFlags update_flags = update_quadrature_points | update_jacobians; - if (spacedim>dim) - update_flags |= update_jacobian_grads; - std_cxx11::unique_ptr mdata (get_data(update_flags, - point_quadrature)); - - mdata->mapping_support_points = this->compute_mapping_support_points (cell); - - // dispatch to the various specializations for spacedim=dim, - // spacedim=dim+1, etc - return do_transform_real_to_unit_cell_internal<1>(cell, p, initial_p_unit, *mdata); -} - -template<> -Point<2> -MappingManifold<2, 2>:: -transform_real_to_unit_cell_internal -(const Triangulation<2, 2>::cell_iterator &cell, - const Point<2> &p, - const Point<2> &initial_p_unit) const -{ - const int dim = 2; - const int spacedim = 2; - - const Quadrature point_quadrature(initial_p_unit); - - UpdateFlags update_flags = update_quadrature_points | update_jacobians; - if (spacedim>dim) - update_flags |= update_jacobian_grads; - std_cxx11::unique_ptr mdata (get_data(update_flags, - point_quadrature)); - - mdata->mapping_support_points = this->compute_mapping_support_points (cell); - - // dispatch to the various specializations for spacedim=dim, - // spacedim=dim+1, etc - return do_transform_real_to_unit_cell_internal<2>(cell, p, initial_p_unit, *mdata); -} - -template<> -Point<3> -MappingManifold<3, 3>:: -transform_real_to_unit_cell_internal -(const Triangulation<3, 3>::cell_iterator &cell, - const Point<3> &p, - const Point<3> &initial_p_unit) const -{ - const int dim = 3; - const int spacedim = 3; - - const Quadrature point_quadrature(initial_p_unit); - - UpdateFlags update_flags = update_quadrature_points | update_jacobians; - if (spacedim>dim) - update_flags |= update_jacobian_grads; - std_cxx11::unique_ptr mdata (get_data(update_flags, - point_quadrature)); - - mdata->mapping_support_points = this->compute_mapping_support_points (cell); - - // dispatch to the various specializations for spacedim=dim, - // spacedim=dim+1, etc - return do_transform_real_to_unit_cell_internal<3>(cell, p, initial_p_unit, *mdata); -} - -template<> -Point<1> -MappingManifold<1, 2>:: -transform_real_to_unit_cell_internal -(const Triangulation<1, 2>::cell_iterator &cell, - const Point<2> &p, - const Point<1> &initial_p_unit) const -{ - const int dim = 1; - const int spacedim = 2; - - const Quadrature point_quadrature(initial_p_unit); - - UpdateFlags update_flags = update_quadrature_points | update_jacobians; - if (spacedim>dim) - update_flags |= update_jacobian_grads; - std_cxx11::unique_ptr mdata (get_data(update_flags, - point_quadrature)); - - mdata->mapping_support_points = this->compute_mapping_support_points (cell); - - // dispatch to the various specializations for spacedim=dim, - // spacedim=dim+1, etc - return do_transform_real_to_unit_cell_internal_codim1<1>(cell, p, initial_p_unit, *mdata); -} - -template<> -Point<2> -MappingManifold<2, 3>:: -transform_real_to_unit_cell_internal -(const Triangulation<2, 3>::cell_iterator &cell, - const Point<3> &p, - const Point<2> &initial_p_unit) const -{ - const int dim = 2; - const int spacedim = 3; - - const Quadrature point_quadrature(initial_p_unit); - - UpdateFlags update_flags = update_quadrature_points | update_jacobians; - if (spacedim>dim) - update_flags |= update_jacobian_grads; - std_cxx11::unique_ptr mdata (get_data(update_flags, - point_quadrature)); - - mdata->mapping_support_points = this->compute_mapping_support_points (cell); - - // dispatch to the various specializations for spacedim=dim, - // spacedim=dim+1, etc - return do_transform_real_to_unit_cell_internal_codim1<2>(cell, p, initial_p_unit, *mdata); -} - -template<> -Point<1> -MappingManifold<1, 3>:: -transform_real_to_unit_cell_internal -(const Triangulation<1, 3>::cell_iterator &, - const Point<3> &, - const Point<1> &) const -{ - Assert (false, ExcNotImplemented()); - return Point<1>(); -} - - - -template -Point -MappingManifold:: -transform_real_to_unit_cell (const typename Triangulation::cell_iterator &cell, - const Point &p) const -{ - // Use an exact formula if one is available. this is only the case - // for Q1 mappings in 1d, and in 2d if dim==spacedim - if ((polynomial_degree == 1) && - ((dim == 1) - || - ((dim == 2) && (dim == spacedim)))) - { - // The dimension-dependent algorithms are much faster (about 25-45x in - // 2D) but fail most of the time when the given point (p) is not in the - // cell. The dimension-independent Newton algorithm given below is - // slower, but more robust (though it still sometimes fails). Therefore - // this function implements the following strategy based on the - // p's dimension: - // - // * In 1D this mapping is linear, so the mapping is always invertible - // (and the exact formula is known) as long as the cell has non-zero - // length. - // * In 2D the exact (quadratic) formula is called first. If either the - // exact formula does not succeed (negative discriminant in the - // quadratic formula) or succeeds but finds a solution outside of the - // unit cell, then the Newton solver is called. The rationale for the - // second choice is that the exact formula may provide two different - // answers when mapping a point outside of the real cell, but the - // Newton solver (if it converges) will only return one answer. - // Otherwise the exact formula successfully found a point in the unit - // cell and that value is returned. - // * In 3D there is no (known to the authors) exact formula, so the Newton - // algorithm is used. - const std_cxx11::array, GeometryInfo::vertices_per_cell> - vertices = this->get_vertices(cell); - try - { - switch (dim) - { - case 1: - { - // formula not subject to any issues in 1d - if (spacedim == 1) - return internal::MappingQ1::transform_real_to_unit_cell(vertices, p); - else - { - const std::vector > a (vertices.begin(), - vertices.end()); - return internal::MappingQ1::transform_real_to_unit_cell_initial_guess(a,p); - } - } - - case 2: - { - const Point point - = internal::MappingQ1::transform_real_to_unit_cell(vertices, p); - - // formula not guaranteed to work for points outside of - // the cell. only take the computed point if it lies - // inside the reference cell - const double eps = 1e-15; - if (-eps <= point(1) && point(1) <= 1 + eps && - -eps <= point(0) && point(0) <= 1 + eps) - { - return point; - } - else - break; - } - - default: - { - // we should get here, based on the if-condition at the top - Assert(false, ExcInternalError()); - } - } - } - catch (const typename Mapping::ExcTransformationFailed &) - { - // simply fall through and continue on to the standard Newton code - } - } - else - { - // we can't use an explicit formula, - } - - - Point initial_p_unit; - if (polynomial_degree == 1) - { - // Find the initial value for the Newton iteration by a normal - // projection to the least square plane determined by the vertices - // of the cell - const std::vector > a - = this->compute_mapping_support_points (cell); - Assert(a.size() == GeometryInfo::vertices_per_cell, - ExcInternalError()); - initial_p_unit = internal::MappingQ1::transform_real_to_unit_cell_initial_guess(a,p); - } - else - { - try - { - // Find the initial value for the Newton iteration by a normal - // projection to the least square plane determined by the vertices - // of the cell - // - // we do this by first getting all support points, then - // throwing away all but the vertices, and finally calling - // the same function as above - std::vector > a - = this->compute_mapping_support_points (cell); - a.resize(GeometryInfo::vertices_per_cell); - initial_p_unit = internal::MappingQ1::transform_real_to_unit_cell_initial_guess(a,p); - } - catch (const typename Mapping::ExcTransformationFailed &) - { - for (unsigned int d=0; d::project_to_unit_cell(initial_p_unit); - } - - // perform the Newton iteration and return the result. note that - // this statement may throw an exception, which we simply pass up to - // the caller - return this->transform_real_to_unit_cell_internal(cell, p, initial_p_unit); -} - - template UpdateFlags @@ -1934,7 +331,7 @@ typename MappingManifold::InternalData * MappingManifold::get_face_data (const UpdateFlags update_flags, const Quadrature &quadrature) const { - InternalData *data = new InternalData(polynomial_degree); + InternalData *data = new InternalData(); data->initialize_face (this->requires_update_flags(update_flags), QProjector::project_to_all_faces(quadrature), quadrature.size()); @@ -1949,7 +346,7 @@ typename MappingManifold::InternalData * MappingManifold::get_subface_data (const UpdateFlags update_flags, const Quadrature& quadrature) const { - InternalData *data = new InternalData(polynomial_degree); + InternalData *data = new InternalData(); data->initialize_face (this->requires_update_flags(update_flags), QProjector::project_to_all_subfaces(quadrature), quadrature.size()); @@ -3641,7 +2038,7 @@ add_line_support_points (const typename Triangulation::cell_iterat cell->get_manifold() : line->get_manifold() ); - get_intermediate_points_on_object (manifold, line_support_points, line, line_points); + // get_intermediate_points_on_object (manifold, line_support_points, line, line_points); if (dim==3) { @@ -3668,134 +2065,134 @@ MappingManifold<3,3>:: add_quad_support_points(const Triangulation<3,3>::cell_iterator &cell, std::vector > &a) const { - const unsigned int faces_per_cell = GeometryInfo<3>::faces_per_cell, - vertices_per_face = GeometryInfo<3>::vertices_per_face, - lines_per_face = GeometryInfo<3>::lines_per_face, - vertices_per_cell = GeometryInfo<3>::vertices_per_cell; - - static const StraightBoundary<3> straight_boundary; - // used if face quad at boundary or entirely in the interior of the domain - std::vector > quad_points ((polynomial_degree-1)*(polynomial_degree-1)); - // used if only one line of face quad is at boundary - std::vector > b(4*polynomial_degree); - - // Used by the new Manifold interface. This vector collects the - // vertices used to compute the intermediate points. - std::vector > vertices(4); - - // loop over all faces and collect points on them - for (unsigned int face_no=0; face_no::face_iterator face = cell->face(face_no); - - // select the correct mappings for the present face - const bool face_orientation = cell->face_orientation(face_no), - face_flip = cell->face_flip (face_no), - face_rotation = cell->face_rotation (face_no); - -#ifdef DEBUG - // some sanity checks up front - for (unsigned int i=0; ivertex_index(i)==cell->vertex_index( - GeometryInfo<3>::face_to_cell_vertices(face_no, i, - face_orientation, - face_flip, - face_rotation)), - ExcInternalError()); - - // indices of the lines that bound a face are given by GeometryInfo<3>:: - // face_to_cell_lines - for (unsigned int i=0; iline(i)==cell->line(GeometryInfo<3>::face_to_cell_lines( - face_no, i, face_orientation, face_flip, face_rotation)), - ExcInternalError()); -#endif - - // if face at boundary, then ask boundary object to return intermediate - // points on it - if (face->at_boundary()) - { - get_intermediate_points_on_object(face->get_manifold(), line_support_points, face, quad_points); - - // in 3D, the orientation, flip and rotation of the face might not - // match what we expect here, namely the standard orientation. thus - // reorder points accordingly. since a Mapping uses the same shape - // function as an FE_Q, we can ask a FE_Q to do the reordering for us. - for (unsigned int i=0; iadjust_quad_dof_index_for_face_orientation(i, - face_orientation, - face_flip, - face_rotation)]); - } - else - { - // face is not at boundary, but maybe some of its lines are. count - // them - unsigned int lines_at_boundary=0; - for (unsigned int i=0; iline(i)->at_boundary()) - ++lines_at_boundary; - - Assert(lines_at_boundary<=lines_per_face, ExcInternalError()); - - // if at least one of the lines bounding this quad is at the - // boundary, then collect points separately - if (lines_at_boundary>0) - { - // call of function add_weighted_interior_points increases size of b - // about 1. There resize b for the case the mentioned function - // was already called. - b.resize(4*polynomial_degree); - - // b is of size 4*degree, make sure that this is the right size - Assert(b.size()==vertices_per_face+lines_per_face*(polynomial_degree-1), - ExcDimensionMismatch(b.size(), - vertices_per_face+lines_per_face*(polynomial_degree-1))); - - // sort the points into b. We used access from the cell (not - // from the face) to fill b, so we can assume a standard face - // orientation. Doing so, the calculated points will be in - // standard orientation as well. - for (unsigned int i=0; i::face_to_cell_vertices(face_no, i)]; - - for (unsigned int i=0; i::face_to_cell_lines( - face_no, i)*(polynomial_degree-1)+j]; - - // Now b includes the support points on the quad and we can - // apply the laplace vector - add_weighted_interior_points (support_point_weights_on_quad, b); - AssertDimension (b.size(), - 4*this->polynomial_degree + - (this->polynomial_degree-1)*(this->polynomial_degree-1)); - - for (unsigned int i=0; i<(polynomial_degree-1)*(polynomial_degree-1); ++i) - a.push_back(b[4*polynomial_degree+i]); - } - else - { - // face is entirely in the interior. get intermediate - // points from the relevant manifold object. - vertices.resize(4); - for (unsigned int i=0; i<4; ++i) - vertices[i] = face->vertex(i); - get_intermediate_points (face->get_manifold(), line_support_points, vertices, quad_points); - // in 3D, the orientation, flip and rotation of the face might - // not match what we expect here, namely the standard - // orientation. thus reorder points accordingly. since a Mapping - // uses the same shape function as an FE_Q, we can ask a FE_Q to - // do the reordering for us. - for (unsigned int i=0; iadjust_quad_dof_index_for_face_orientation(i, - face_orientation, - face_flip, - face_rotation)]); - } - } - } +// const unsigned int faces_per_cell = GeometryInfo<3>::faces_per_cell, +// vertices_per_face = GeometryInfo<3>::vertices_per_face, +// lines_per_face = GeometryInfo<3>::lines_per_face, +// vertices_per_cell = GeometryInfo<3>::vertices_per_cell; + +// static const StraightBoundary<3> straight_boundary; +// // used if face quad at boundary or entirely in the interior of the domain +// std::vector > quad_points ((polynomial_degree-1)*(polynomial_degree-1)); +// // used if only one line of face quad is at boundary +// std::vector > b(4*polynomial_degree); + +// // Used by the new Manifold interface. This vector collects the +// // vertices used to compute the intermediate points. +// std::vector > vertices(4); + +// // loop over all faces and collect points on them +// for (unsigned int face_no=0; face_no::face_iterator face = cell->face(face_no); + +// // select the correct mappings for the present face +// const bool face_orientation = cell->face_orientation(face_no), +// face_flip = cell->face_flip (face_no), +// face_rotation = cell->face_rotation (face_no); + +// #ifdef DEBUG +// // some sanity checks up front +// for (unsigned int i=0; ivertex_index(i)==cell->vertex_index( +// GeometryInfo<3>::face_to_cell_vertices(face_no, i, +// face_orientation, +// face_flip, +// face_rotation)), +// ExcInternalError()); + +// // indices of the lines that bound a face are given by GeometryInfo<3>:: +// // face_to_cell_lines +// for (unsigned int i=0; iline(i)==cell->line(GeometryInfo<3>::face_to_cell_lines( +// face_no, i, face_orientation, face_flip, face_rotation)), +// ExcInternalError()); +// #endif + +// // if face at boundary, then ask boundary object to return intermediate +// // points on it +// if (face->at_boundary()) +// { +// get_intermediate_points_on_object(face->get_manifold(), line_support_points, face, quad_points); + +// // in 3D, the orientation, flip and rotation of the face might not +// // match what we expect here, namely the standard orientation. thus +// // reorder points accordingly. since a Mapping uses the same shape +// // function as an FE_Q, we can ask a FE_Q to do the reordering for us. +// for (unsigned int i=0; iadjust_quad_dof_index_for_face_orientation(i, +// face_orientation, +// face_flip, +// face_rotation)]); +// } +// else +// { +// // face is not at boundary, but maybe some of its lines are. count +// // them +// unsigned int lines_at_boundary=0; +// for (unsigned int i=0; iline(i)->at_boundary()) +// ++lines_at_boundary; + +// Assert(lines_at_boundary<=lines_per_face, ExcInternalError()); + +// // if at least one of the lines bounding this quad is at the +// // boundary, then collect points separately +// if (lines_at_boundary>0) +// { +// // call of function add_weighted_interior_points increases size of b +// // about 1. There resize b for the case the mentioned function +// // was already called. +// b.resize(4*polynomial_degree); + +// // b is of size 4*degree, make sure that this is the right size +// Assert(b.size()==vertices_per_face+lines_per_face*(polynomial_degree-1), +// ExcDimensionMismatch(b.size(), +// vertices_per_face+lines_per_face*(polynomial_degree-1))); + +// // sort the points into b. We used access from the cell (not +// // from the face) to fill b, so we can assume a standard face +// // orientation. Doing so, the calculated points will be in +// // standard orientation as well. +// for (unsigned int i=0; i::face_to_cell_vertices(face_no, i)]; + +// for (unsigned int i=0; i::face_to_cell_lines( +// face_no, i)*(polynomial_degree-1)+j]; + +// // Now b includes the support points on the quad and we can +// // apply the laplace vector +// add_weighted_interior_points (support_point_weights_on_quad, b); +// AssertDimension (b.size(), +// 4*this->polynomial_degree + +// (this->polynomial_degree-1)*(this->polynomial_degree-1)); + +// for (unsigned int i=0; i<(polynomial_degree-1)*(polynomial_degree-1); ++i) +// a.push_back(b[4*polynomial_degree+i]); +// } +// else +// { +// // face is entirely in the interior. get intermediate +// // points from the relevant manifold object. +// vertices.resize(4); +// for (unsigned int i=0; i<4; ++i) +// vertices[i] = face->vertex(i); +// get_intermediate_points (face->get_manifold(), line_support_points, vertices, quad_points); +// // in 3D, the orientation, flip and rotation of the face might +// // not match what we expect here, namely the standard +// // orientation. thus reorder points accordingly. since a Mapping +// // uses the same shape function as an FE_Q, we can ask a FE_Q to +// // do the reordering for us. +// for (unsigned int i=0; iadjust_quad_dof_index_for_face_orientation(i, +// face_orientation, +// face_flip, +// face_rotation)]); +// } +// } +// } } @@ -3806,11 +2203,11 @@ MappingManifold<2,3>:: add_quad_support_points(const Triangulation<2,3>::cell_iterator &cell, std::vector > &a) const { - std::vector > quad_points ((polynomial_degree-1)*(polynomial_degree-1)); - get_intermediate_points_on_object (cell->get_manifold(), line_support_points, - cell, quad_points); - for (unsigned int i=0; i > quad_points ((polynomial_degree-1)*(polynomial_degree-1)); + // get_intermediate_points_on_object (cell->get_manifold(), line_support_points, + // cell, quad_points); + // for (unsigned int i=0; i