From: Marco Tezzele Date: Thu, 2 Apr 2015 12:48:28 +0000 (+0200) Subject: new wonderful maping X-Git-Tag: v8.3.0-rc1~316^2~1 X-Git-Url: https://gitweb.dealii.org/cgi-bin/gitweb.cgi?a=commitdiff_plain;h=ea9efd738b9e82e60585018328dfcd3471def25a;p=dealii.git new wonderful maping --- diff --git a/include/deal.II/fe/mapping_fe.h b/include/deal.II/fe/mapping_fe.h new file mode 100644 index 0000000000..a5e21c1058 --- /dev/null +++ b/include/deal.II/fe/mapping_fe.h @@ -0,0 +1,742 @@ +// --------------------------------------------------------------------- +// $Id: mapping_fe.h 30450 2013-08-23 15:48:29Z kronbichler $ +// +// Copyright (C) 2001 - 2013 by the deal.II authors +// +// This file is part of the deal.II library. +// +// The deal.II library is free software; you can use it, redistribute +// it, and/or modify it under the terms of the GNU Lesser General +// Public License as published by the Free Software Foundation; either +// version 2.1 of the License, or (at your option) any later version. +// The full text of the license can be found in the file LICENSE at +// the top level of the deal.II distribution. +// +// --------------------------------------------------------------------- + +#ifndef __deal2__mapping_fe_h +#define __deal2__mapping_fe_h + + +#include +#include +#include +#include +#include +#include + + +DEAL_II_NAMESPACE_OPEN + + +/*!@addtogroup mapping */ +/*@{*/ + +/** + * The MappingFE is a generalization of the MappingQEulerian class, for arbitrary + * vectorial finite elements. The main difference is that this class uses a vector + * of absolute positions instead of a vector of displacement. + * In particular we think of a collections of a FE_Q or + * Bezier finite element (FE_Bernstein) repeated a number of times equal to the space + * dimension. The idea is to construct the mapping using a vector of control + * points, a DoFHandler associated to the geometry of the problem and a + * ComponentMask that tells us which components to use for the mapping. + * This mapping will grab from the DoFHandler the finite element, or better + * the collection of finite elements, to compute the mapping shape functions. + * So we will have two different Finite Element and DoFHandler, one for the + * solution field and one to describe the geometry of the problem. Historically + * in the deal.II library there was not such a distinction. The differences + * between this mapping and the MappingQ class are quite important. + * The MappingFE, being a generalization, requires a higher level of abstraction. + * This is the reason why it takes a DoFHandler and a vector of control points + * that are the coefficients of the shape function (so in general it is a vector + * of coefficient). + * + * + * Typically, the DoFHandler operates on a finite element that is constructed + * as a system element (FESystem) from continuous FE_Q() objects. An example + * is shown below: + * @code + * const FE_Q feq(1); + * const FESystem fesystem(feq, spacedim); + * DoFHandler dhq(triangulation); + * dhq.distribute_dofs(fesystem); + * Vector eulerq(dhq.n_dofs()); + * const ComponentMask mask(spacedim, true); + * MappingFE map(eulerq, dhq, mask); + * map.update_euler_vector_using_triangulation(eulerq); + * @endcode + + + * + * @author Luca Heltai, Marco Tezzele 2013, 2015 + */ +template , + class VECTOR=Vector > +class MappingFE : public Mapping +{ +public: + /** + * Constructor. The first argument is a VECTOR that specifies the + * transformation of the domain from the reference to the current + * configuration. This is filled calling the method + * update_euler_vector_using_triangulation. + */ + MappingFE (const VECTOR &euler_vector, + const DH &euler_dof_handler, + const ComponentMask mask=ComponentMask()); + + /** + * Copy constructor. Performs a deep copy, i.e. duplicates what #tensor_pols + * points to instead of simply copying the #tensor_pols pointer as done by a + * default copy constructor. + */ + MappingFE (const MappingFE &mapping); + + /** + * Destructor. + */ + virtual ~MappingFE (); + + /** Fill the euler vector with + the information coming from + the triangulation. Makes this + map equivalent to MappingQ1, + and it works ONLY if the + underlying fe has support + points. */ + void update_euler_vector_using_triangulation(VECTOR &vector); + + + + /** + * Transforms the point @p p on the unit cell to the point @p p_real on the + * real cell @p cell and returns @p p_real. + */ + virtual Point + transform_unit_to_real_cell ( + const typename Triangulation::cell_iterator &cell, + const Point &p) const; + + /** + * Transforms the point @p p on the real cell to the point @p p_unit on the + * unit cell @p cell and returns @p p_unit. + * + * Uses Newton iteration and the @p transform_unit_to_real_cell function. + * + * In the codimension one case, this function returns the normal projection + * of the real point @p p on the curve or surface identified by the @p cell. + */ + virtual Point + transform_real_to_unit_cell ( + const typename Triangulation::cell_iterator &cell, + const Point &p) const; + + + virtual void + transform (const VectorSlice > > input, + VectorSlice > > output, + const typename Mapping::InternalDataBase &internal, + const MappingType type) const; + + virtual void + transform (const VectorSlice > > input, + VectorSlice > > output, + const typename Mapping::InternalDataBase &internal, + const MappingType type) const; + + virtual + void + transform (const VectorSlice > > input, + VectorSlice > > output, + const typename Mapping::InternalDataBase &internal, + const MappingType type) const; + + + + /** + * Return the degree of the mapping, i.e. the value which was passed to the + * constructor. + */ + unsigned int get_degree () const; + + /** + * Return the ComponentMask of the mapping, i.e. which components to use for + * the mapping. + */ + ComponentMask get_fe_mask () const; + + /** + * Return a pointer to a copy of the present object. The caller of this copy + * then assumes ownership of it. + */ + virtual + Mapping *clone () const; + + + /** + * Storage for internal data of + * d-linear transformation. + */ + class InternalData : public Mapping::InternalDataBase + { + public: + /** + * Constructor. + */ + InternalData(const FiniteElement &fe, + const ComponentMask mask); + + /** + * Shape function at quadrature + * point. Shape functions are + * in tensor product order, so + * vertices must be reordered + * to obtain transformation. + */ + 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. + */ + 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. + */ + 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); + + /** + * Return an estimate (in + * bytes) or the memory + * consumption of this + * object. + */ + 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. + * + * 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; + + /** + * 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. + */ + 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. + */ + std::vector< DerivativeForm<1,dim,spacedim> > contravariant; + + /** + * 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; + + /** + * Auxiliary vectors for internal use. + */ + std::vector > > aux; + + /** + * 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. + */ + unsigned int n_shape_functions; + + ComponentMask mask; + }; + + + /** + * Transforms a point @p p on + * the unit cell to the point + * @p p_real on the real cell + * @p cell and returns @p p_real. + * + * This function is called by + * @p transform_unit_to_real_cell + * and multiple times (through the + * Newton iteration) by + * @p transform_real_to_unit_cell_internal. + * + * Takes a reference to an + * @p InternalData that must + * already include the shape + * values at point @p p and the + * mapping support points of the + * cell. + * + * This @p InternalData argument + * avoids multiple computations + * of the shape values at point + * @p p and especially multiple + * computations of the mapping + * support points. + */ + Point + transform_unit_to_real_cell_internal (const InternalData &mdata) const; + + + /** + * Transforms the point @p p on + * the real cell to the corresponding + * point on the unit cell + * @p cell by a Newton + * iteration. + * + * Takes a reference to an + * @p InternalData that is + * assumed to be previously + * created by the @p get_data + * function with @p UpdateFlags + * including + * @p update_transformation_values + * and + * @p update_transformation_gradients + * and a one point Quadrature + * that includes the given + * initial guess for the + * transformation + * @p initial_p_unit. Hence this + * function assumes that + * @p mdata already includes the + * transformation shape values + * and gradients computed at + * @p initial_p_unit. + * + * @p mdata will be changed by + * this function. + */ + Point + transform_real_to_unit_cell_internal (const typename Triangulation::cell_iterator &cell, + const Point &p, + const Point &initial_p_unit, + InternalData &mdata) const; + + /** + * Do the computation for the + * fill_* functions. + */ + void compute_fill (const typename Triangulation::cell_iterator &cell, + const unsigned int npts, + const typename QProjector::DataSetDescriptor data_set, + const CellSimilarity::Similarity cell_similarity, + InternalData &data, + std::vector > &quadrature_points) const; + + + /** + * Do the computation for the + * fill_* functions. + */ + void compute_fill_face (const typename Triangulation::cell_iterator &cell, + const unsigned int face_no, + const unsigned int subface_no, + const unsigned int npts, + const typename QProjector::DataSetDescriptor data_set, + const std::vector &weights, + InternalData &mapping_data, + std::vector > &quadrature_points, + std::vector &JxW_values, + std::vector > &boundary_form, + std::vector > &normal_vectors, + std::vector > &jacobians, + std::vector > &inverse_jacobians) const; + + + /** + * Always returns @p false. + */ + virtual + bool preserves_vertex_locations () const; + + DeclException0(ExcInactiveCell); + +protected: + /** + * Implementation of the interface in Mapping. + */ + virtual void + fill_fe_values (const typename Triangulation::cell_iterator &cell, + const Quadrature &quadrature, + typename Mapping::InternalDataBase &mapping_data, + typename std::vector > &quadrature_points, + std::vector &JxW_values, + std::vector > &jacobians, + std::vector > &jacobian_grads, + std::vector > &inverse_jacobians, + std::vector > &cell_normal_vectors, + CellSimilarity::Similarity &cell_similarity) const ; + + /** + * Implementation of the interface in Mapping. + */ + virtual void + fill_fe_face_values (const typename Triangulation::cell_iterator &cell, + const unsigned int face_no, + const Quadrature& quadrature, + typename Mapping::InternalDataBase &mapping_data, + typename std::vector > &quadrature_points, + std::vector &JxW_values, + std::vector > &exterior_forms, + std::vector > &normal_vectors, + std::vector > &jacobians, + std::vector > &inverse_jacobians) const ; + + /** + * Implementation of the interface in Mapping. + */ + virtual void + fill_fe_subface_values (const typename Triangulation::cell_iterator &cell, + const unsigned int face_no, + const unsigned int sub_no, + const Quadrature& quadrature, + typename Mapping::InternalDataBase &mapping_data, + typename std::vector > &quadrature_points, + std::vector &JxW_values, + std::vector > &exterior_forms, + std::vector > &normal_vectors, + std::vector > &jacobians, + std::vector > &inverse_jacobians) const ; + + + /** + This function and the next allow to generate the transform require by + the virtual transform() in mapping, but unfortunately in C++ one cannot + declare a virtual template function. + */ + template < int rank > + void + transform_fields(const VectorSlice > > input, + VectorSlice< std::vector > > output, + const typename Mapping::InternalDataBase &internal, + const MappingType type) const; + + + /** + see doc in transform_fields + */ + template < int rank > + void + transform_differential_forms( + const VectorSlice > > input, + VectorSlice > > output, + const typename Mapping::InternalDataBase &mapping_data, + const MappingType mapping_type) const; + + +protected: + /** + * Reference to the vector of shifts. + */ + + SmartPointer >euler_vector; + /** + * A FiniteElement 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. We could make this a pointer to prevent + * construction in 1D and 2D, but since memory and time requirements are not + * particularly high this seems unnecessary at the moment. + */ + SmartPointer, MappingFE > fe; + + + /** + * Pointer to the DoFHandler to which the mapping vector is associated. + */ + SmartPointer >euler_dof_handler; + + + +private: +// + /** + * Update internal degrees of + * freedom. */ + void update_internal_dofs(const typename Triangulation::cell_iterator &cell) const; + + + mutable std::vector local_dofs; + + mutable std::vector dof_indices; + + /** + * Mutex to protect local_dofs. + */ + + mutable Threads::Mutex mutex; + + + virtual void + compute_shapes_virtual (const std::vector > &unit_points, + typename MappingFE::InternalData &data) const; + + UpdateFlags + update_once (const UpdateFlags in) const; + + UpdateFlags + update_each (const UpdateFlags in) const; + + void + compute_data (const UpdateFlags update_flags, + const Quadrature &q, + const unsigned int n_original_q_points, + InternalData &data) const; + + void + compute_face_data (const UpdateFlags update_flags, + const Quadrature &q, + const unsigned int n_original_q_points, + InternalData &data) const; + + virtual + typename Mapping::InternalDataBase * + get_data (const UpdateFlags, + const Quadrature &quadrature) const; + + virtual + typename Mapping::InternalDataBase * + get_face_data (const UpdateFlags flags, + const Quadrature& quadrature) const; + + virtual + typename Mapping::InternalDataBase * + get_subface_data (const UpdateFlags flags, + const Quadrature& quadrature) const; + + + /* + * Which components to use for the mapping. + */ + const ComponentMask fe_mask; + + + /** + * Mapping between indices in the FE space and the real space. This vector contains one + * index for each component of the finite element space. If the index is one for which + * the ComponentMask which is used to construct this element is false, then + * numbers::invalid_unsigned_int is returned, otherwise the component in real space is + * returned. For example, if we construct the mapping using ComponentMask(spacedim, true), + * then this vector contains {0,1,2} in spacedim = 3. + */ + std::vector fe_to_real; + + + /** + * Declare other MappingFE classes friends. + */ + template friend class MappingFE; +}; + +/*@}*/ + +/* -------------- declaration of explicit specializations ------------- */ + +#ifndef DOXYGEN + + +template +inline +double +MappingFE::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 & +MappingFE::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 +Tensor<1,dim> +MappingFE::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> & +MappingFE::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 +Tensor<2,dim> +MappingFE::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> & +MappingFE::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 +bool +MappingFE::preserves_vertex_locations () const +{ + return false; +} + + + + +#endif // DOXYGEN + +DEAL_II_NAMESPACE_CLOSE + +#endif diff --git a/source/fe/CMakeLists.txt b/source/fe/CMakeLists.txt index b676e48baf..80cfbcaa60 100644 --- a/source/fe/CMakeLists.txt +++ b/source/fe/CMakeLists.txt @@ -48,6 +48,7 @@ SET(_src mapping_c1.cc mapping_cartesian.cc mapping.cc + mapping_fe.cc mapping_q1.cc mapping_q1_eulerian.cc mapping_q.cc @@ -87,6 +88,7 @@ SET(_inst mapping_c1.inst.in mapping_cartesian.inst.in mapping.inst.in + mapping_fe.inst.in mapping_q1_eulerian.inst.in mapping_q1.inst.in mapping_q_eulerian.inst.in diff --git a/source/fe/mapping_fe.cc b/source/fe/mapping_fe.cc new file mode 100644 index 0000000000..e66dd6ec7a --- /dev/null +++ b/source/fe/mapping_fe.cc @@ -0,0 +1,1386 @@ +// --------------------------------------------------------------------- +// $Id: mapping_fe.cc 30450 2013-08-23 15:48:29Z kronbichler $ +// +// Copyright (C) 2001 - 2013 by the deal.II authors +// +// This file is part of the deal.II library. +// +// The deal.II library is free software; you can use it, redistribute +// it, and/or modify it under the terms of the GNU Lesser General +// Public License as published by the Free Software Foundation; either +// version 2.1 of the License, or (at your option) any later version. +// The full text of the license can be found in the file LICENSE at +// the top level of the deal.II distribution. +// +// --------------------------------------------------------------------- + +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include + +#include +#include +#include + + + +DEAL_II_NAMESPACE_OPEN + + +template +MappingFE::InternalData::InternalData (const FiniteElement &fe, + const ComponentMask mask) + : + n_shape_functions (fe.dofs_per_cell), + mask (mask) +{} + + + +template +std::size_t +MappingFE::InternalData::memory_consumption () const +{ + return 0; +} + + +template +MappingFE::MappingFE (const VECTOR &euler_vector, + const DH &euler_dof_handler, + const ComponentMask mask) + : + euler_vector(&euler_vector), + fe(&euler_dof_handler.get_fe()), + euler_dof_handler(&euler_dof_handler), + local_dofs(fe->dofs_per_cell), + dof_indices(fe->dofs_per_cell), + fe_mask(mask.size() ? mask : + ComponentMask(fe->get_nonzero_components(0).size(), true)), + fe_to_real(fe_mask.size(), numbers::invalid_unsigned_int) +{ + + unsigned int size = 0; + for (unsigned int i=0; i +MappingFE::MappingFE (const MappingFE &mapping) + : + euler_vector(mapping.euler_vector), + fe(mapping.fe), + euler_dof_handler(mapping.euler_dof_handler), + local_dofs(mapping.local_dofs), + dof_indices(mapping.dof_indices), + fe_mask(mapping.fe_mask), + fe_to_real(mapping.fe_to_real) +{} + + +template +MappingFE::~MappingFE () +{ + euler_dof_handler = NULL; + fe = NULL; + euler_vector = NULL; +} + + +template +void +MappingFE::compute_shapes_virtual ( + const std::vector > &unit_points, + typename MappingFE::InternalData &data) const +{ + const unsigned int n_points=unit_points.size(); + + if (data.shape_values.size()!=0 || data.shape_derivatives.size()!=0) + for (unsigned int point=0; pointshape_value(i, unit_points[point]); + + if (data.shape_derivatives.size()!=0) + for (unsigned int i=0; ishape_grad(i, unit_points[point]); + + if (data.shape_second_derivatives.size()!=0) + for (unsigned int i=0; ishape_grad_grad(i, unit_points[point]); + } +} + + +template +UpdateFlags +MappingFE::update_once (const UpdateFlags in) const +{ + UpdateFlags out = UpdateFlags(in & (update_transformation_values + | update_transformation_gradients)); + + // Shape function values + if (in & update_quadrature_points) + out |= update_transformation_values; + + // Shape function gradients + if (in & (update_covariant_transformation + | update_contravariant_transformation + | update_JxW_values + | update_boundary_forms + | update_normal_vectors + | update_jacobians + | update_jacobian_grads + | update_inverse_jacobians)) + out |= update_transformation_gradients; + + return out; +} + + + +template +UpdateFlags +MappingFE::update_each (const UpdateFlags in) const +{ + // Select flags of concern for the + // transformation. + UpdateFlags out = UpdateFlags(in & (update_quadrature_points + | update_covariant_transformation + | update_contravariant_transformation + | update_JxW_values + | update_boundary_forms + | update_normal_vectors + | update_volume_elements + | update_jacobians + | update_jacobian_grads + | update_inverse_jacobians)); + + // add flags if the respective + // quantities are necessary to + // compute what we need. note that + // some flags appear in both + // conditions and in subsequents + // set operations. this leads to + // some circular logic. the only + // way to treat this is to + // iterate. since there are 4 + // if-clauses in the loop, it will + // take at most 3 iterations to + // converge. do them: + for (unsigned int i=0; i<4; ++i) + { + // The following is a little incorrect: + // If not applied on a face, + // update_boundary_forms does not + // make sense. On the other hand, + // it is necessary on a + // face. Currently, + // update_boundary_forms is simply + // ignored for the interior of a + // cell. + if (out & (update_JxW_values + | update_normal_vectors)) + out |= update_boundary_forms; + + if (out & (update_covariant_transformation + | update_JxW_values + | update_jacobians + | update_jacobian_grads + | update_boundary_forms + | update_normal_vectors)) + out |= update_contravariant_transformation; + + if (out & (update_inverse_jacobians)) + out |= update_covariant_transformation; + + // The contravariant transformation + // is a Piola transformation, which + // requires the determinant of the + // Jacobi matrix of the transformation. + // Therefore these values have to be + // updated for each cell. + if (out & update_contravariant_transformation) + out |= update_JxW_values; + + if (out & update_normal_vectors) + out |= update_JxW_values; + } + + return out; +} + + + + +template +void +MappingFE::compute_data (const UpdateFlags update_flags, + const Quadrature &q, + const unsigned int n_original_q_points, + InternalData &data) const +{ + const unsigned int n_q_points = q.size(); + + data.update_once = update_once(update_flags); + data.update_each = update_each(update_flags); + data.update_flags = data.update_once | data.update_each; + + const UpdateFlags flags(data.update_flags); + + if (flags & update_transformation_values) + data.shape_values.resize(data.n_shape_functions * n_q_points); + + if (flags & update_transformation_gradients) + data.shape_derivatives.resize(data.n_shape_functions * n_q_points); + + if (flags & update_covariant_transformation) + data.covariant.resize(n_original_q_points); + + if (flags & update_contravariant_transformation) + data.contravariant.resize(n_original_q_points); + + if (flags & update_volume_elements) + data.volume_elements.resize(n_original_q_points); + + if (flags & update_jacobian_grads) + data.shape_second_derivatives.resize(data.n_shape_functions * n_q_points); + + compute_shapes_virtual (q.get_points(), data); + +} + + +template +void +MappingFE::compute_face_data (const UpdateFlags update_flags, + const Quadrature &q, + const unsigned int n_original_q_points, + InternalData &data) const +{ + compute_data (update_flags, q, n_original_q_points, data); + + if (dim > 1) + { + if (data.update_flags & update_boundary_forms) + { + data.aux.resize (dim-1, std::vector > (n_original_q_points)); + + // Compute tangentials to the + // unit cell. + const unsigned int nfaces = GeometryInfo::faces_per_cell; + data.unit_tangentials.resize (nfaces*(dim-1), + std::vector > (n_original_q_points)); + if (dim==2) + { + // ensure a counterclock wise + // 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 (data.unit_tangentials[i].begin(), + data.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 (data.unit_tangentials[i].begin(), + data.unit_tangentials[i].end(), tang1); + std::fill (data.unit_tangentials[nfaces+i].begin(), + data.unit_tangentials[nfaces+i].end(), tang2); + } + } + } + } +} + + +template +typename Mapping::InternalDataBase * +MappingFE::get_data (const UpdateFlags update_flags, + const Quadrature &quadrature) const +{ + InternalData *data = new InternalData(*fe, fe_mask); + this->compute_data (update_flags, quadrature, + quadrature.size(), *data); + return data; +} + + + +template +typename Mapping::InternalDataBase * +MappingFE::get_face_data (const UpdateFlags update_flags, + const Quadrature& quadrature) const +{ + InternalData *data = new InternalData(*fe, fe_mask); + const Quadrature q (QProjector::project_to_all_faces(quadrature)); + this->compute_face_data (update_flags, q, + quadrature.size(), *data); + + return data; +} + + +template +typename Mapping::InternalDataBase * +MappingFE::get_subface_data (const UpdateFlags update_flags, + const Quadrature& quadrature) const +{ + InternalData *data = new InternalData(*fe, fe_mask); + const Quadrature q (QProjector::project_to_all_subfaces(quadrature)); + this->compute_face_data (update_flags, q, + quadrature.size(), *data); + + return data; +} + + +// Note that the CellSimilarity flag is modifyable, since MappingFE can need to +// recalculate data even when cells are similar. +template +void +MappingFE::fill_fe_values ( + const typename Triangulation::cell_iterator &cell, + const Quadrature &q, + typename Mapping::InternalDataBase &mapping_data, + std::vector > &quadrature_points, + std::vector &JxW_values, + std::vector > &jacobians, + std::vector > &jacobian_grads, + std::vector > &inverse_jacobians, + std::vector > &normal_vectors, + CellSimilarity::Similarity &cell_similarity) const +{ + AssertDimension(fe->dofs_per_cell, local_dofs.size()); + Assert(local_dofs.size()>0, ExcMessage("Cannot do anything with zero degrees of freedom")); + + // convert data object to internal data for this class. fails with an + // exception if that is not possible + Assert (dynamic_cast (&mapping_data) != 0, ExcInternalError()); + InternalData &data = static_cast (mapping_data); + + // depending on this result, use this or the other data object for the + // mapping. furthermore, we need to ensure that the flag indicating whether + // we can use some similarity has to be modified - for a general MappingFE, + // the data needs to be recomputed anyway since then the mapping changes the + // data. this needs to be known also for later operations, so modify the + // variable here. this also affects the calculation of the next cell -- if + // we use Q1 data on the next cell, the data will still be invalid. + + if (get_degree() > 1) + cell_similarity = CellSimilarity::invalid_next_cell; + + const unsigned int n_q_points=q.size(); + + compute_fill (cell, n_q_points, QProjector::DataSetDescriptor::cell (), cell_similarity, + data, quadrature_points); + + const UpdateFlags update_flags(data.current_update_flags()); + const std::vector &weights=q.get_weights(); + + // Multiply quadrature weights by absolute value of Jacobian determinants or + // the area element g=sqrt(DX^t DX) in case of codim > 0 + + if (update_flags & (update_normal_vectors | update_JxW_values)) + { + AssertDimension (JxW_values.size(), n_q_points); + + Assert( !(update_flags & update_normal_vectors ) || + (normal_vectors.size() == n_q_points), + ExcDimensionMismatch(normal_vectors.size(), n_q_points)); + + + if (cell_similarity != CellSimilarity::translation) + for (unsigned int point=0; point 1e-12*Utilities::fixed_power(cell->diameter()/ + std::sqrt(double(dim))), + (typename Mapping::ExcDistortedMappedCell(cell->center(), det, point))); + JxW_values[point] = weights[point] * det; + } + // if dim==spacedim, then there is no cell normal to + // compute. since this is for FEValues (and not FEFaceValues), + // there are also no face normals to compute + else //codim>0 case + { + Tensor<1, spacedim> DX_t [dim]; + for (unsigned int i=0; i G; //First fundamental form + for (unsigned int i=0; idirection_flag() == false) + normal_vectors[point] *= -1.; + } + + } + } //codim>0 case + + } + } + + + // copy values from InternalData to vector given by reference + if (update_flags & update_jacobians) + { + AssertDimension (jacobians.size(), n_q_points); + if (cell_similarity != CellSimilarity::translation) + for (unsigned int point=0; point()); + + const unsigned int data_set = QProjector::DataSetDescriptor::cell(); + + for (unsigned int point=0; point *second = + &data.second_derivative(point+data_set, 0); + + double result [spacedim][dim][dim]; + + for (unsigned int k=0; ksystem_to_component_index(k).first; + if (fe_mask[comp_k]) + for (unsigned int j=0; j +void +MappingFE::fill_fe_face_values ( + const typename Triangulation::cell_iterator &cell, + const unsigned int face_no, + const Quadrature &q, + typename Mapping::InternalDataBase &mapping_data, + std::vector > &quadrature_points, + std::vector &JxW_values, + std::vector > &exterior_forms, + std::vector > &normal_vectors, + std::vector > &jacobians, + std::vector > &inverse_jacobians) const +// std::vector > &exterior_forms, +// std::vector > &normal_vectors) const +{ + // convert data object to internal data for this class. fails with an + // exception if that is not possible + +// AssertThrow(false, ExcNotImplemented()); + + + + Assert (dynamic_cast (&mapping_data) != 0, + ExcInternalError()); + InternalData &data = static_cast (mapping_data); + + const unsigned int n_q_points=q.size(); + this->compute_fill_face (cell, face_no, numbers::invalid_unsigned_int, + n_q_points, + QProjector::DataSetDescriptor:: + face (face_no, + cell->face_orientation(face_no), + cell->face_flip(face_no), + cell->face_rotation(face_no), + n_q_points), + q.get_weights(), + data, + quadrature_points, JxW_values, + exterior_forms, normal_vectors, jacobians, + inverse_jacobians); + // quadrature_points, JxW_values, + // exterior_forms, normal_vectors); +} + + +template +void +MappingFE::fill_fe_subface_values (const typename Triangulation::cell_iterator &cell, + const unsigned int face_no, + const unsigned int sub_no, + const Quadrature &q, + typename Mapping::InternalDataBase &mapping_data, + std::vector > &quadrature_points, + std::vector &JxW_values, + std::vector > &exterior_forms, + std::vector > &normal_vectors, + std::vector > &jacobians, + std::vector > &inverse_jacobians) const +// std::vector > &exterior_forms, +// std::vector > &normal_vectors) const +{ + //AssertThrow(false, ExcNotImplemented()); + + + // convert data object to internal data for this class. fails with an + // exception if that is not possible + Assert (dynamic_cast (&mapping_data) != 0, + ExcInternalError()); + InternalData &data = static_cast (mapping_data); + + const unsigned int n_q_points=q.size(); + this->compute_fill_face (cell, face_no, sub_no, + n_q_points, + QProjector::DataSetDescriptor:: + subface (face_no, sub_no, + cell->face_orientation(face_no), + cell->face_flip(face_no), + cell->face_rotation(face_no), + n_q_points, + cell->subface_case(face_no)), + q.get_weights(), + data, + quadrature_points, JxW_values, + exterior_forms, normal_vectors, jacobians, + inverse_jacobians); + // quadrature_points, JxW_values, + // exterior_forms, normal_vectors); +} + + + + +template +void +MappingFE::transform ( + const VectorSlice > > input, + VectorSlice > > output, + const typename Mapping::InternalDataBase &mapping_data, + const MappingType mapping_type) const +{ + AssertDimension (input.size(), output.size()); + + transform_fields(input, output, mapping_data, mapping_type); +} + + + +template +void +MappingFE::transform ( + const VectorSlice > > input, + VectorSlice > > output, + const typename Mapping::InternalDataBase &mapping_data, + const MappingType mapping_type) const +{ + AssertDimension (input.size(), output.size()); + + std::vector > aux_output1(output.size()); + VectorSlice< std::vector > > aux_output( aux_output1); + + transform_differential_forms(input, aux_output, mapping_data, mapping_type); + + for (unsigned int i=0; i +void MappingFE::transform +(const VectorSlice > > input, + VectorSlice > > output, + const typename Mapping::InternalDataBase &mapping_data, + const MappingType mapping_type) const +{ + AssertDimension (input.size(), output.size()); + + AssertThrow(false, ExcNotImplemented()); +} + + +template +template < int rank > +void MappingFE::transform_fields( + const VectorSlice > > input, + VectorSlice > > output, + const typename Mapping::InternalDataBase &mapping_data, + const MappingType mapping_type) const +{ + AssertDimension (input.size(), output.size()); + Assert (dynamic_cast(&mapping_data) != 0, + ExcInternalError()); + const InternalData &data = static_cast(mapping_data); + + switch (mapping_type) + { + case mapping_contravariant: + { + Assert (data.update_flags & update_contravariant_transformation, + typename FEValuesBase::ExcAccessToUninitializedField("update_contravariant_transformation")); + + for (unsigned int i=0; i::ExcAccessToUninitializedField("update_contravariant_transformation")); + Assert (data.update_flags & update_volume_elements, + typename FEValuesBase::ExcAccessToUninitializedField("update_volume_elements")); + Assert (rank==1, ExcMessage("Only for rank 1")); + for (unsigned int i=0; i::ExcAccessToUninitializedField("update_contravariant_transformation")); + + for (unsigned int i=0; i +template < int rank > +void MappingFE::transform_differential_forms( + const VectorSlice > > input, + VectorSlice > > output, + const typename Mapping::InternalDataBase &mapping_data, + const MappingType mapping_type) const +{ + + AssertDimension (input.size(), output.size()); + Assert (dynamic_cast(&mapping_data) != 0, + ExcInternalError()); + const InternalData &data = static_cast(mapping_data); + + switch (mapping_type) + { + case mapping_covariant: + { + Assert (data.update_flags & update_contravariant_transformation, + typename FEValuesBase::ExcAccessToUninitializedField("update_contravariant_transformation")); + + for (unsigned int i=0; i +Point +MappingFE:: +transform_unit_to_real_cell (const typename Triangulation::cell_iterator &cell, + const Point &p) const +{ +// Use the get_data function to create an InternalData with data vectors of +// the right size and transformation shape values already computed at point +// p. Make sure only one processor at a time performs this operations. + + Threads::Mutex::ScopedLock lock(mutex); + update_internal_dofs(cell); + const Quadrature point_quadrature(p); + std::auto_ptr + mdata (dynamic_cast ( + get_data(update_transformation_values, point_quadrature))); + + return this->transform_unit_to_real_cell_internal(*mdata); +} + + +template +Point +MappingFE:: +transform_unit_to_real_cell_internal (const InternalData &data) const +{ + Point p_real; + + for (unsigned int i=0; isystem_to_component_index(i).first; + if (fe_mask[comp_i]) + p_real[fe_to_real[comp_i]] += local_dofs[i] * data.shape(0,i); + } + + return p_real; +} + + + +template +Point +MappingFE:: +transform_real_to_unit_cell (const typename Triangulation::cell_iterator &cell, + const Point &p) const +{ + { + Threads::Mutex::ScopedLock lock(mutex); + update_internal_dofs(cell); + } + // first a Newton iteration based on the real mapping. It uses the center + // point of the cell as a starting point + Point initial_p_unit; + + + for (unsigned int d=0; d::project_to_unit_cell(initial_p_unit); + + const Quadrature point_quadrature(initial_p_unit); + + UpdateFlags update_flags = update_transformation_values|update_transformation_gradients; + if (spacedim>dim) + update_flags |= update_jacobian_grads; + std::auto_ptr + mdata (dynamic_cast ( + get_data(update_flags,point_quadrature))); + + return this->transform_real_to_unit_cell_internal(cell, p, initial_p_unit, *mdata); + +} + + +template +Point +MappingFE:: +transform_real_to_unit_cell_internal +(const typename Triangulation::cell_iterator &cell, + const Point &p, + const Point &initial_p_unit, + InternalData &mdata) const +{ + { + Threads::Mutex::ScopedLock lock(mutex); + + update_internal_dofs(cell); + } + + const unsigned int n_shapes=mdata.shape_values.size(); + Assert(n_shapes!=0, ExcInternalError()); + AssertDimension (mdata.shape_derivatives.size(), n_shapes); + + + // Newton iteration to solve + // f(x)=p(x)-p=0 + // x_{n+1}=x_n-[f'(x)]^{-1}f(x) + + // The start value was set to be the + // linear approximation to the cell + + // The shape values and derivatives + // of the mapping at this point are + // previously computed. + + // For the <2,3> case there is a + // template specialization. + + // f(x) + + //Point p_minus_F; + Point p_unit = initial_p_unit; + Point f; + + compute_shapes_virtual(std::vector > (1, p_unit), mdata); + Point p_real(transform_unit_to_real_cell_internal(mdata)); + + Tensor<1,spacedim> p_minus_F = p - p_real; + + const double eps = 1.e-12*cell->diameter(); + const unsigned int newton_iteration_limit = 20; + + unsigned int newton_iteration=0; + + while (p_minus_F.norm_square() > eps*eps) + { + // f'(x) + Point DF[dim]; + Tensor<2,dim> df; + + for (unsigned int k=0; k &grad_k = mdata.derivative(0,k); + + unsigned int comp_k = fe->system_to_component_index(k).first; + if (fe_mask[comp_k]) + for (unsigned int j=0; j delta; + contract (delta, invert(df), static_cast&>(f)); + + // do a line search + double step_length = 1; + do + { + // update of p_unit. The + // spacedimth 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), mdata); + + // f(x) + Point p_real_trial = transform_unit_to_real_cell_internal(mdata); + //const Point f_trial = p - p_real_trial; + const Tensor<1,spacedim> f_trial = p - p_real_trial; + + // see if we are making progress with the current step length + // and if not, reduce it by a factor of two and try again + if (f_trial.norm() < p_minus_F.norm()) + { + p_real = p_real_trial; + p_unit = p_unit_trial; + + p_minus_F = f_trial; + break; + } + else if (step_length > 0.05) + step_length /= 2; + else + { + std::cout << "Line search failed. With dim = " << dim << " spacedim = " + << spacedim << std::endl; + goto failure; + } + } + while (true); + + ++newton_iteration; + if (newton_iteration > newton_iteration_limit) + { + std::cout << "Too many newton iterations. With dim = " << dim << " spacedim = " + << spacedim << std::endl; + goto failure; + } + } + + return p_unit; + + // if we get to the following label, then we have either run out + // of Newton iterations, or the line search has not converged. + // in either case, we need to give up, so throw an exception that + // can then be caught +failure: + AssertThrow (false, (typename Mapping::ExcTransformationFailed())); + + // ...the compiler wants us to return something, though we can + // of course never get here... + return Point(); +} + + +template +void +MappingFE::compute_fill ( + const typename Triangulation::cell_iterator &cell, + const unsigned int n_q_points, + const typename QProjector::DataSetDescriptor data_set, + const CellSimilarity::Similarity cell_similarity, + InternalData &data, + std::vector > &quadrature_points) const +{ + const UpdateFlags update_flags(data.current_update_flags()); + Threads::Mutex::ScopedLock lock(mutex); + update_internal_dofs(cell); + + // first compute quadrature points + if (update_flags & update_quadrature_points) + { + AssertDimension (quadrature_points.size(), n_q_points); + + for (unsigned int point=0; point result; + const double *shape = &data.shape(point+data_set,0); + + for (unsigned int k=0; ksystem_to_component_index(k).first; + if (fe_mask[comp_k]) + result[fe_to_real[comp_k]] += local_dofs[k] * shape[k]; + } + + quadrature_points[point] = result; + } + } + + + // then Jacobians + if (update_flags & update_contravariant_transformation) + { + AssertDimension (data.contravariant.size(), n_q_points); + + // if the current cell is just a + // translation of the previous one, no + // need to recompute jacobians... + if (cell_similarity != CellSimilarity::translation) + { + std::fill(data.contravariant.begin(), data.contravariant.end(), + DerivativeForm<1,dim,spacedim>()); + + Assert (data.n_shape_functions > 0, ExcInternalError()); + + for (unsigned int point=0; point *data_derv = + &data.derivative(point+data_set, 0); + + Tensor<1, dim> result[spacedim]; + + for (unsigned int k=0; ksystem_to_component_index(k).first; + if (fe_mask[comp_k]) + result[fe_to_real[comp_k]] += local_dofs[k] * data_derv[k]; + } + + // write result into contravariant data. for + // j=dim in the case dim +void +MappingFE::compute_fill_face ( + const typename Triangulation::cell_iterator &cell, + const unsigned int face_no, + const unsigned int subface_no, + const unsigned int n_q_points,//npts + const typename QProjector::DataSetDescriptor data_set, + const std::vector &weights, + InternalData &data, + std::vector > &quadrature_points, + std::vector &JxW_values, + std::vector > &boundary_forms, + std::vector > &normal_vectors, + std::vector > &jacobians, + std::vector > &inverse_jacobians) const +{ + compute_fill (cell, n_q_points, data_set, CellSimilarity::none, + data, quadrature_points); + + + const UpdateFlags update_flags(data.current_update_flags()); + + if (update_flags & update_boundary_forms) + { + AssertDimension (boundary_forms.size(), n_q_points); + if (update_flags & update_normal_vectors) + AssertDimension (normal_vectors.size(), n_q_points); + if (update_flags & update_JxW_values) + AssertDimension (JxW_values.size(), n_q_points); + + // map the unit tangentials to the real cell. checking for d!=dim-1 + // eliminates compiler warnings regarding unsigned int expressions < + // 0. + for (unsigned int d=0; d!=dim-1; ++d) + { + Assert (face_no+GeometryInfo::faces_per_cell*d < + data.unit_tangentials.size(), + ExcInternalError()); + Assert (data.aux[d].size() <= + data.unit_tangentials[face_no+GeometryInfo::faces_per_cell*d].size(), + ExcInternalError()); + + transform (data.unit_tangentials[face_no+GeometryInfo::faces_per_cell*d], + data.aux[d], + data, + mapping_contravariant); + } + + // if dim==spacedim, we can use the unit tangentials to compute the + // boundary form by simply taking the cross product + if (dim == spacedim) + { + for (unsigned int i=0; i cell_normal; + const DerivativeForm<1,spacedim,dim> DX_t = + data.contravariant[point].transpose(); + cross_product(cell_normal,DX_t[0],DX_t[1]); + cell_normal /= cell_normal.norm(); + + // then compute the face normal from the face tangent + // and the cell normal: + cross_product (boundary_forms[point], + data.aux[0][point], cell_normal); + + } + + } + } + + + + if (update_flags & (update_normal_vectors + | update_JxW_values)) + for (unsigned int i=0; i::subface_ratio( + cell->subface_case(face_no), subface_no); + JxW_values[i] *= area_ratio; + } + } + + if (update_flags & update_normal_vectors) + normal_vectors[i] = Point(boundary_forms[i] / boundary_forms[i].norm()); + } + } + +} + + +template +unsigned int +MappingFE::get_degree() const +{ + return fe->degree; +} + + +template +ComponentMask +MappingFE::get_fe_mask() const +{ + return this->fe_mask; +} + + +template +Mapping * +MappingFE::clone () const +{ + return new MappingFE(*this); +} + + +template +void +MappingFE::update_internal_dofs ( + const typename Triangulation::cell_iterator &cell) const +{ + if (euler_dof_handler == 0) + { + std::cout << "euler_dof_handler is empty!" << std::endl; + return; + } + + typename DH::cell_iterator dof_cell(*cell, euler_dof_handler); + Assert (dof_cell->active() == true, ExcInactiveCell()); + + dof_cell->get_dof_indices(dof_indices); + + for (unsigned int i=0; i +void MappingFE::update_euler_vector_using_triangulation +(VECTOR &vector) +{ + if ( fe->has_support_points() ) + { + std::vector > support_points = fe->get_unit_support_points(); + typename DH::active_cell_iterator cell; + Quadrature quad(support_points); + + MappingQ map_q(fe->degree); + FEValues fe_v(map_q, *fe, quad, update_quadrature_points); + std::vector dofs(fe->dofs_per_cell); + + AssertDimension(fe->dofs_per_cell, support_points.size()); + Assert(fe->is_primitive(), ExcMessage("FE is not Primitive! This won't work.")); + + for (cell = euler_dof_handler->begin_active(); cell != euler_dof_handler->end(); ++cell) + { + fe_v.reinit(cell); + cell->get_dof_indices(dofs); + const std::vector > &points = fe_v.get_quadrature_points(); + for (unsigned int q = 0; q < points.size(); ++q) + { + unsigned int comp = fe->system_to_component_index(q).first; + vector(dofs[q]) = points[q][comp]; + } + } + + } + else + { + // Construct a MappingFE with an FEQ + FESystem feq(FE_Q(fe->degree), spacedim); + DH dhq(euler_dof_handler->get_tria()); + dhq.distribute_dofs(feq); + VECTOR eulerq(dhq.n_dofs()); + const ComponentMask maskq(spacedim, true); + MappingFE newfe(eulerq, dhq, maskq); + + newfe.update_euler_vector_using_triangulation(eulerq); + + FullMatrix transfer(fe->dofs_per_cell, feq.dofs_per_cell); + std::vector > points = feq.get_unit_support_points(); + + // Here construct the matrix!!!! + for (unsigned int i=0; idofs_per_cell; ++i) + { + for (unsigned int j=0; jsystem_to_component_index(i).first + == + feq.system_to_component_index(j).first) + transfer(j,i) = fe->shape_value(i, points[j]); + } + } + VectorTools::interpolate(dhq, *euler_dof_handler, transfer, eulerq, vector); + } +} + + + + +// explicit instantiations +#include "mapping_fe.inst" + + +DEAL_II_NAMESPACE_CLOSE diff --git a/source/fe/mapping_fe.inst.in b/source/fe/mapping_fe.inst.in new file mode 100644 index 0000000000..a06af31506 --- /dev/null +++ b/source/fe/mapping_fe.inst.in @@ -0,0 +1,26 @@ +// --------------------------------------------------------------------- +// $Id: mapping_fe.inst.in 30049 2013-07-18 19:42:40Z maier $ +// +// Copyright (C) 1998 - 2013 by the deal.II authors +// +// This file is part of the deal.II library. +// +// The deal.II library is free software; you can use it, redistribute +// it, and/or modify it under the terms of the GNU Lesser General +// Public License as published by the Free Software Foundation; either +// version 2.1 of the License, or (at your option) any later version. +// The full text of the license can be found in the file LICENSE at +// the top level of the deal.II distribution. +// +// --------------------------------------------------------------------- + + + + +for (deal_II_dimension : DIMENSIONS; deal_II_space_dimension : SPACE_DIMENSIONS) + { +#if deal_II_dimension <= deal_II_space_dimension + template class MappingFE, dealii::Vector >; +#endif + } + diff --git a/tests/fe/mapping_fe_real_to_unit_q1.cc b/tests/fe/mapping_fe_real_to_unit_q1.cc new file mode 100644 index 0000000000..10097bb9f5 --- /dev/null +++ b/tests/fe/mapping_fe_real_to_unit_q1.cc @@ -0,0 +1,141 @@ +// --------------------------------------------------------------------- +// $Id$ +// +// Copyright (C) 2006 - 2013 by the deal.II authors +// +// This file is part of the deal.II library. +// +// The deal.II library is free software; you can use it, redistribute +// it, and/or modify it under the terms of the GNU Lesser General +// Public License as published by the Free Software Foundation; either +// version 2.1 of the License, or (at your option) any later version. +// The full text of the license can be found in the file LICENSE at +// the top level of the deal.II distribution. +// +// --------------------------------------------------------------------- + + + +// on a somewhat deformed cube, verify that if we push forward a bunch +// of points from the reference to the real cell and then call +// MappingFE::transform_unit_to_real_cell that we get the same point as +// we had in the beginning. + +#include "../tests.h" + +#include +#include +#include +#include +#include +#include +#include +#include + +using namespace dealii; + + +template +void test_real_to_unit_cell() +{ + deallog << "dim=" << dim << ", spacedim=" << spacedim << std::endl; + + Triangulation triangulation; + GridGenerator::hyper_cube (triangulation); + + const unsigned int n_points = 4; + std::vector< Point > unit_points(Utilities::fixed_power(n_points)); + + switch (dim) + { + case 1: + for (unsigned int x=0; x feq(1); + const FESystem fesystem(feq, spacedim); + DoFHandler dhq(triangulation); + dhq.distribute_dofs(fesystem); + Vector eulerq(dhq.n_dofs()); + const ComponentMask mask(spacedim, true); + + MappingFE map(eulerq, dhq, mask); + + map.update_euler_vector_using_triangulation(eulerq); + + typename Triangulation::active_cell_iterator + cell = triangulation.begin_active(); + + //Move a vertex a little bit + const unsigned int n_dx = 5; + const double dx = 0.4/n_dx; + Point direction; + for (unsigned int j=0; jvertex(0) = double(j)*direction; + + for (unsigned int i=0; i p = map.transform_unit_to_real_cell(cell,unit_points[i]); + const Point p_unit = map.transform_real_to_unit_cell(cell,p); + + Assert (unit_points[i].distance(p_unit) < 1e-10, ExcInternalError()); + } + } + + deallog << "OK" << std::endl; +} + + +int +main() +{ + std::ofstream logfile ("output"); + deallog.attach(logfile); + deallog.depth_console(0); + deallog.threshold_double(1.e-10); + + test_real_to_unit_cell<1,1>(); + test_real_to_unit_cell<2,2>(); + test_real_to_unit_cell<3,3>(); + + test_real_to_unit_cell<1,2>(); + test_real_to_unit_cell<1,3>(); + test_real_to_unit_cell<2,3>(); + + return 0; +} diff --git a/tests/fe/mapping_fe_real_to_unit_q1.output b/tests/fe/mapping_fe_real_to_unit_q1.output new file mode 100644 index 0000000000..c68c1fb460 --- /dev/null +++ b/tests/fe/mapping_fe_real_to_unit_q1.output @@ -0,0 +1,43 @@ + +DEAL::dim=1, spacedim=1 +DEAL::Vertex displacement: 0.00000 +DEAL::Vertex displacement: 0.0800000 +DEAL::Vertex displacement: 0.160000 +DEAL::Vertex displacement: 0.240000 +DEAL::Vertex displacement: 0.320000 +DEAL::OK +DEAL::dim=2, spacedim=2 +DEAL::Vertex displacement: 0.00000 0.00000 +DEAL::Vertex displacement: 0.0800000 0.0800000 +DEAL::Vertex displacement: 0.160000 0.160000 +DEAL::Vertex displacement: 0.240000 0.240000 +DEAL::Vertex displacement: 0.320000 0.320000 +DEAL::OK +DEAL::dim=3, spacedim=3 +DEAL::Vertex displacement: 0.00000 0.00000 0.00000 +DEAL::Vertex displacement: 0.0800000 0.0800000 0.0800000 +DEAL::Vertex displacement: 0.160000 0.160000 0.160000 +DEAL::Vertex displacement: 0.240000 0.240000 0.240000 +DEAL::Vertex displacement: 0.320000 0.320000 0.320000 +DEAL::OK +DEAL::dim=1, spacedim=2 +DEAL::Vertex displacement: 0.00000 0.00000 +DEAL::Vertex displacement: 0.0800000 0.0800000 +DEAL::Vertex displacement: 0.160000 0.160000 +DEAL::Vertex displacement: 0.240000 0.240000 +DEAL::Vertex displacement: 0.320000 0.320000 +DEAL::OK +DEAL::dim=1, spacedim=3 +DEAL::Vertex displacement: 0.00000 0.00000 0.00000 +DEAL::Vertex displacement: 0.0800000 0.0800000 0.0800000 +DEAL::Vertex displacement: 0.160000 0.160000 0.160000 +DEAL::Vertex displacement: 0.240000 0.240000 0.240000 +DEAL::Vertex displacement: 0.320000 0.320000 0.320000 +DEAL::OK +DEAL::dim=2, spacedim=3 +DEAL::Vertex displacement: 0.00000 0.00000 0.00000 +DEAL::Vertex displacement: 0.0800000 0.0800000 0.0800000 +DEAL::Vertex displacement: 0.160000 0.160000 0.160000 +DEAL::Vertex displacement: 0.240000 0.240000 0.240000 +DEAL::Vertex displacement: 0.320000 0.320000 0.320000 +DEAL::OK diff --git a/tests/fe/mapping_fe_real_to_unit_q5_curved.cc b/tests/fe/mapping_fe_real_to_unit_q5_curved.cc new file mode 100644 index 0000000000..486c306991 --- /dev/null +++ b/tests/fe/mapping_fe_real_to_unit_q5_curved.cc @@ -0,0 +1,151 @@ +// --------------------------------------------------------------------- +// $Id$ +// +// Copyright (C) 2006 - 2013 by the deal.II authors +// +// This file is part of the deal.II library. +// +// The deal.II library is free software; you can use it, redistribute +// it, and/or modify it under the terms of the GNU Lesser General +// Public License as published by the Free Software Foundation; either +// version 2.1 of the License, or (at your option) any later version. +// The full text of the license can be found in the file LICENSE at +// the top level of the deal.II distribution. +// +// --------------------------------------------------------------------- + + + +// on a somewhat deformed cube, verify that if we push forward a bunch +// of points from the reference to the real cell and then call +// MappingFE::transform_unit_to_real_cell that we get the same point as +// we had in the beginning. + +// We use a Q5 mapping but this time we +// actually curve one boundary of the cell which ensures that the +// mapping is really higher order than just Q1 + +#include "../tests.h" + +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include + +#include +#include + +using namespace dealii; + + +template +void test_real_to_unit_cell() +{ + deallog << "dim=" << dim << ", spacedim=" << spacedim << std::endl; + + // define a boundary that fits the + // the vertices of the hyper cube + // we're going to create below + HyperBallBoundary boundary (Point(), + std::sqrt(1.*dim)); + + Triangulation triangulation; + GridGenerator::hyper_cube (triangulation, -1, 1); + + // set the boundary indicator for + // one face of the single cell + triangulation.set_boundary (1, boundary); + triangulation.begin_active()->face(0)->set_boundary_indicator (1); + + + const unsigned int n_points = 5; + std::vector< Point > unit_points(Utilities::fixed_power(n_points)); + + switch (dim) + { + case 1: + for (unsigned int x=0; x feq(5); + const FESystem fesystem(feq, spacedim); + DoFHandler dhq(triangulation); + dhq.distribute_dofs(fesystem); + Vector eulerq(dhq.n_dofs()); + const ComponentMask mask(spacedim, true); + + MappingFE map(eulerq, dhq, mask); + + map.update_euler_vector_using_triangulation(eulerq); + + + typename Triangulation::active_cell_iterator + cell = triangulation.begin_active(); + + for (unsigned int i=0; i p = map.transform_unit_to_real_cell(cell,unit_points[i]); + + const Point p_unit = map.transform_real_to_unit_cell(cell,p); + + AssertThrow (unit_points[i].distance(p_unit) < 1e-10, ExcInternalError()); + } + + deallog << "OK" << std::endl; + +} + + +int +main() +{ + std::ofstream logfile ("output"); + deallog.attach(logfile); + deallog.depth_console(0); + deallog.threshold_double(1.e-10); + + test_real_to_unit_cell<1,1>(); + test_real_to_unit_cell<2,2>(); + test_real_to_unit_cell<3,3>(); + + test_real_to_unit_cell<1,2>(); + test_real_to_unit_cell<2,3>(); + + + //test_real_to_unit_cell<1,3>(); + return 0; +} diff --git a/tests/fe/mapping_fe_real_to_unit_q5_curved.output b/tests/fe/mapping_fe_real_to_unit_q5_curved.output new file mode 100644 index 0000000000..5693e8cff9 --- /dev/null +++ b/tests/fe/mapping_fe_real_to_unit_q5_curved.output @@ -0,0 +1,11 @@ + +DEAL::dim=1, spacedim=1 +DEAL::OK +DEAL::dim=2, spacedim=2 +DEAL::OK +DEAL::dim=3, spacedim=3 +DEAL::OK +DEAL::dim=1, spacedim=2 +DEAL::OK +DEAL::dim=2, spacedim=3 +DEAL::OK