From: Luca Heltai Date: Mon, 6 Apr 2015 18:06:32 +0000 (+0200) Subject: Renamed MappingFE to MappingFEField X-Git-Tag: v8.3.0-rc1~267^2~4 X-Git-Url: https://gitweb.dealii.org/cgi-bin/gitweb.cgi?a=commitdiff_plain;h=22b6d1e6367700d6591ef974402b6be85bb498dd;p=dealii.git Renamed MappingFE to MappingFEField Fixed documentation. Fixed indentation. Inverted DH and VECTOR in template declaration. Renamed get_fe_mask to get_component_mask Added entry in changes.h fixed: destructor, useless lines, commented code. Renamed tests added missing documentation and solved the thread safe issue eliminated if statement in compute_shapes_virtual solved all the problems highlighted in https://github.com/dealii/dealii/pull/727 const quadrature points in update_euler_vector_using_triangulation Added two test with a non trivial component mask (mapping_fe_field_real_to_unit_q1_mask.cc) and fixed if statement in update_euler_vector_using_triangulation Added assert dimension in update_euler_vector_using_triangulation start working on else statement in update_euler added one test with fe_bernstein (still failing) Added tests from q1 to q5 and from b1 to b5 Fixed doxygen comments. Added a few comments and fixed a const reference. fixed inital guess in transform_real_to_unit added hyperballboundary in test b2_mask --- diff --git a/doc/news/changes.h b/doc/news/changes.h index aae3579f89..2c3d4c119b 100644 --- a/doc/news/changes.h +++ b/doc/news/changes.h @@ -511,6 +511,12 @@ inconvenience this causes. (Wolfgang Bangerth, 2015/04/08) +
  • New: A new MappingFEField() class has been added to the library + that generalizes MappingQEulerian to allow arbitrary FiniteElements. +
    + (Marco Tezzele, Luca Heltai, 2015/04/06) +
  • +
  • Changed: In the spirit of the changes made to the distinction between Point and Tensor objects discussed above, the first argument to GridTools::shift() has been changed from a Point to a Tensor@<1,dim@>. diff --git a/include/deal.II/fe/mapping_fe.h b/include/deal.II/fe/mapping_fe_field.h similarity index 72% rename from include/deal.II/fe/mapping_fe.h rename to include/deal.II/fe/mapping_fe_field.h index 308c0a82db..7e065af27c 100644 --- a/include/deal.II/fe/mapping_fe.h +++ b/include/deal.II/fe/mapping_fe_field.h @@ -32,29 +32,29 @@ DEAL_II_NAMESPACE_OPEN /*@{*/ /** - * 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 displacements. - * 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). + * The MappingFEField is a generalization of the MappingQEulerian + * class, for arbitrary vector finite elements. The two main + * differences are that this class uses a vector of absolute positions + * instead of a vector of displacements, and it allows for arbitrary + * FiniteElement types, instead of only FE_Q. * + * This class effectively decouples the topology from the geometry, by + * relegating all geometrical information to some components of a + * FiniteElement vector field. The components that are used for the + * geometry can be arbitrarily selected at construction time. + * + * The idea is to consider the Triangulation as a parameter + * configuration space, on which we construct an arbitrary + * geometrical mapping, using the instruments of the deal.II library: + * a vector of degrees of freedom, a DoFHandler associated to the + * geometry of the problem and a ComponentMask that tells us which + * components of the FiniteElement to use for the mapping. + * + * Typically, the DoFHandler operates on a finite element that is + * constructed as a system element (FESystem) from continuous FE_Q() + * (for iso-parametric discretizations) or FE_Bernstein() (for + * iso-geometric discretizations) objects. An example is shown below: * - * 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); @@ -62,50 +62,86 @@ DEAL_II_NAMESPACE_OPEN * dhq.distribute_dofs(fesystem); * Vector eulerq(dhq.n_dofs()); * const ComponentMask mask(spacedim, true); - * MappingFE map(eulerq, dhq, mask); + * MappingFEField 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 + class VECTOR=Vector, + class DH=DoFHandler > +class MappingFEField : 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. + * configuration. + * + * In general this class decouples geometry from topology, allowing + * users to define geometries which are only topologically + * equivalent to the underlying Triangulation, but which may + * otherwise be arbitrary. Differently from what happens in + * MappingQEulerian, the FiniteElement field which is passed to the + * constructor is interpreted as an absolute geometrical + * configuration, therefore one has to make sure that the + * euler_vector actually represents a valid geometry (i.e., one with + * no inverted cells, or with no zero-volume cells). + * + * In order to facilitate the user, we provide an + * update_euler_vector_using_triangulation() method to initialize a + * euler_vector in a way similar to what happens in + * MappingQEulerian, by creating a euler_vector which interpolates + * the underlying Triangulation. + * + * If the underlying FiniteElement is a system of FE_Q(), and + * euler_vector is initialized using + * update_euler_vector_using_triangulation(), then this class is in + * all respects identical to MappingQ(). + * + * The optional ComponentMask argument can be used to specify what + * components of the FiniteElement to use for the geometrical + * transformation. If no mask is specified at construction time, + * then a default one is used, which makes this class works in the + * same way of MappingQEulerian(), i.e., the first spacedim + * components of the FiniteElement are assumed to represent the + * geometry of the problem. + * + * Notice that if a mask is specified, it has to match in size the + * underlying FiniteElement, and it has to have exactly spacedim + * non-zero elements, indicating the components (in order) of the + * FiniteElement which will be used for the geometry. + * + * If an incompatible mask is passed, an exception is thrown. */ - MappingFE (const VECTOR &euler_vector, - const DH &euler_dof_handler, - const ComponentMask mask=ComponentMask()); + MappingFEField (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); + MappingFEField (const MappingFEField &mapping); /** - * Destructor. + * Helper function to fill the given euler vector with information + * coming from the triangulation, by interpolating the geometry of + * the Triangulation associated with the DoFHandler used at + * construction time. + * + * The resulting map is guaranteed to be interpolatory at the + * vertices of the Triangulation. Notice that this may or may not be + * meaningful, depending on the FiniteElement you have used to + * construct this MappingFEField. + * + * If the underlying FiniteElement is a system of FE_Q(), and + * euler_vector is initialized using this function, then this class + * is in all respects identical to MappingQ(). */ - 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); + void update_euler_vector_using_triangulation(VECTOR &vector) const; @@ -132,19 +168,27 @@ public: const typename Triangulation::cell_iterator &cell, const Point &p) const; - + /** Reimplemented from Mapping. See the documentation of the base class for + * detailed information. + */ virtual void transform (const VectorSlice > > input, VectorSlice > > output, const typename Mapping::InternalDataBase &internal, const MappingType type) const; + /** Reimplemented from Mapping. See the documentation of the base class for + * detailed information. + */ virtual void transform (const VectorSlice > > input, VectorSlice > > output, const typename Mapping::InternalDataBase &internal, const MappingType type) const; + /** Reimplemented from Mapping. See the documentation of the base class for + * detailed information. + */ virtual void transform (const VectorSlice > > input, @@ -164,7 +208,7 @@ public: * Return the ComponentMask of the mapping, i.e. which components to use for * the mapping. */ - ComponentMask get_fe_mask () const; + ComponentMask get_component_mask () const; /** * Return a pointer to a copy of the present object. The caller of this copy @@ -336,20 +380,26 @@ public: 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. + * 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; + /** + * Stores the mask given at construction time. If no mask was + * specified at construction time, then a default one is used, + * which makes this class works in the same way of + * MappingQEulerian(), i.e., the first spacedim components of the + * FiniteElement are used for the euler_vector and the euler_dh. + * + * If a mask is specified, then it has to match the underlying + * FiniteElement, and it has to have exactly spacedim non-zero + * elements, indicating the components (in order) of the + * FiniteElement which will be used for the euler vector and the + * euler dof handler. + */ ComponentMask mask; }; @@ -537,7 +587,7 @@ protected: * Reference to the vector of shifts. */ - SmartPointer >euler_vector; + 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 @@ -545,13 +595,13 @@ protected: * construction in 1D and 2D, but since memory and time requirements are not * particularly high this seems unnecessary at the moment. */ - SmartPointer, MappingFE > fe; + SmartPointer, MappingFEField > fe; /** * Pointer to the DoFHandler to which the mapping vector is associated. */ - SmartPointer >euler_dof_handler; + SmartPointer >euler_dof_handler; @@ -563,49 +613,78 @@ private: void update_internal_dofs(const typename Triangulation::cell_iterator &cell) const; + /** Reimplemented from Mapping. See the documentation of the base class for + * detailed information. + */ mutable std::vector local_dofs; + /** Reimplemented from Mapping. See the documentation of the base class for + * detailed information. + */ mutable std::vector dof_indices; /** * Mutex to protect local_dofs. */ - mutable Threads::Mutex mutex; + /** Reimplemented from Mapping. See the documentation of the base class for + * detailed information. + */ virtual void compute_shapes_virtual (const std::vector > &unit_points, - typename MappingFE::InternalData &data) const; + typename MappingFEField::InternalData &data) const; + /** Reimplemented from Mapping. See the documentation of the base class for + * detailed information. + */ UpdateFlags update_once (const UpdateFlags in) const; + /** Reimplemented from Mapping. See the documentation of the base class for + * detailed information. + */ UpdateFlags update_each (const UpdateFlags in) const; + /** Reimplemented from Mapping. See the documentation of the base class for + * detailed information. + */ void compute_data (const UpdateFlags update_flags, const Quadrature &q, const unsigned int n_original_q_points, InternalData &data) const; + /** Reimplemented from Mapping. See the documentation of the base class for + * detailed information. + */ void compute_face_data (const UpdateFlags update_flags, const Quadrature &q, const unsigned int n_original_q_points, InternalData &data) const; + /** Reimplemented from Mapping. See the documentation of the base class for + * detailed information. + */ virtual typename Mapping::InternalDataBase * get_data (const UpdateFlags, const Quadrature &quadrature) const; + /** Reimplemented from Mapping. See the documentation of the base class for + * detailed information. + */ virtual typename Mapping::InternalDataBase * get_face_data (const UpdateFlags flags, const Quadrature& quadrature) const; + /** Reimplemented from Mapping. See the documentation of the base class for + * detailed information. + */ virtual typename Mapping::InternalDataBase * get_subface_data (const UpdateFlags flags, @@ -630,9 +709,9 @@ private: /** - * Declare other MappingFE classes friends. + * Declare other MappingFEField classes friends. */ - template friend class MappingFE; + template friend class MappingFEField; }; /*@}*/ @@ -645,8 +724,8 @@ private: template inline double -MappingFE::InternalData::shape (const unsigned int qpoint, - const unsigned int shape_nr) const +MappingFEField::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, @@ -659,8 +738,8 @@ MappingFE::InternalData::shape (const unsigned int qpoin template inline double & -MappingFE::InternalData::shape (const unsigned int qpoint, - const unsigned int shape_nr) +MappingFEField::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, @@ -672,7 +751,7 @@ MappingFE::InternalData::shape (const unsigned int qpoin template inline Tensor<1,dim> -MappingFE::InternalData::derivative (const unsigned int qpoint, +MappingFEField::InternalData::derivative (const unsigned int qpoint, const unsigned int shape_nr) const { Assert(qpoint*n_shape_functions + shape_nr < shape_derivatives.size(), @@ -686,7 +765,7 @@ MappingFE::InternalData::derivative (const unsigned int template inline Tensor<1,dim> & -MappingFE::InternalData::derivative (const unsigned int qpoint, +MappingFEField::InternalData::derivative (const unsigned int qpoint, const unsigned int shape_nr) { Assert(qpoint*n_shape_functions + shape_nr < shape_derivatives.size(), @@ -699,7 +778,7 @@ MappingFE::InternalData::derivative (const unsigned int template inline Tensor<2,dim> -MappingFE::InternalData::second_derivative (const unsigned int qpoint, +MappingFEField::InternalData::second_derivative (const unsigned int qpoint, const unsigned int shape_nr) const { Assert(qpoint*n_shape_functions + shape_nr < shape_second_derivatives.size(), @@ -713,7 +792,7 @@ MappingFE::InternalData::second_derivative (const unsign template inline Tensor<2,dim> & -MappingFE::InternalData::second_derivative (const unsigned int qpoint, +MappingFEField::InternalData::second_derivative (const unsigned int qpoint, const unsigned int shape_nr) { Assert(qpoint*n_shape_functions + shape_nr < shape_second_derivatives.size(), @@ -726,7 +805,7 @@ MappingFE::InternalData::second_derivative (const unsign template inline bool -MappingFE::preserves_vertex_locations () const +MappingFEField::preserves_vertex_locations () const { return false; } diff --git a/source/fe/CMakeLists.txt b/source/fe/CMakeLists.txt index def0a6685a..9d7c756970 100644 --- a/source/fe/CMakeLists.txt +++ b/source/fe/CMakeLists.txt @@ -49,7 +49,7 @@ SET(_src mapping_c1.cc mapping_cartesian.cc mapping.cc - mapping_fe.cc + mapping_fe_field.cc mapping_q1.cc mapping_q1_eulerian.cc mapping_q.cc @@ -90,7 +90,7 @@ SET(_inst mapping_c1.inst.in mapping_cartesian.inst.in mapping.inst.in - mapping_fe.inst.in + mapping_fe_field.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_field.cc similarity index 80% rename from source/fe/mapping_fe.cc rename to source/fe/mapping_fe_field.cc index c8343e5f03..c2a199fc73 100644 --- a/source/fe/mapping_fe.cc +++ b/source/fe/mapping_fe_field.cc @@ -27,7 +27,7 @@ #include #include #include -#include +#include #include #include #include @@ -44,8 +44,8 @@ DEAL_II_NAMESPACE_OPEN -template -MappingFE::InternalData::InternalData (const FiniteElement &fe, +template +MappingFEField::InternalData::InternalData (const FiniteElement &fe, const ComponentMask mask) : n_shape_functions (fe.dofs_per_cell), @@ -54,18 +54,18 @@ MappingFE::InternalData::InternalData (const FiniteEleme -template +template std::size_t -MappingFE::InternalData::memory_consumption () const +MappingFEField::InternalData::memory_consumption () const { return 0; } -template -MappingFE::MappingFE (const VECTOR &euler_vector, - const DH &euler_dof_handler, - const ComponentMask mask) +template +MappingFEField::MappingFEField (const VECTOR &euler_vector, + const DH &euler_dof_handler, + const ComponentMask mask) : euler_vector(&euler_vector), fe(&euler_dof_handler.get_fe()), @@ -76,7 +76,6 @@ MappingFE::MappingFE (const VECTOR &euler_vector, 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 (const VECTOR &euler_vector, } -template -MappingFE::MappingFE (const MappingFE &mapping) +template +MappingFEField::MappingFEField (const MappingFEField &mapping) : euler_vector(mapping.euler_vector), fe(mapping.fe), @@ -100,44 +99,34 @@ MappingFE::MappingFE (const MappingFE -MappingFE::~MappingFE () -{ - euler_dof_handler = NULL; - fe = NULL; - euler_vector = NULL; -} - - -template +template void -MappingFE::compute_shapes_virtual ( +MappingFEField::compute_shapes_virtual ( const std::vector > &unit_points, - typename MappingFE::InternalData &data) const + typename MappingFEField::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]); - } + 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 +template UpdateFlags -MappingFE::update_once (const UpdateFlags in) const +MappingFEField::update_once (const UpdateFlags in) const { UpdateFlags out = UpdateFlags(in & (update_transformation_values | update_transformation_gradients)); @@ -162,9 +151,9 @@ MappingFE::update_once (const UpdateFlags in) const -template +template UpdateFlags -MappingFE::update_each (const UpdateFlags in) const +MappingFEField::update_each (const UpdateFlags in) const { // Select flags of concern for the // transformation. @@ -187,11 +176,11 @@ MappingFE::update_each (const UpdateFlags in) const // set operations. this leads to // some circular logic. the only // way to treat this is to - // iterate. since there are 4 + // iterate. since there are 5 // if-clauses in the loop, it will - // take at most 3 iterations to + // take at most 4 iterations to // converge. do them: - for (unsigned int i=0; i<4; ++i) + for (unsigned int i=0; i<5; ++i) { // The following is a little incorrect: // If not applied on a face, @@ -236,12 +225,12 @@ MappingFE::update_each (const UpdateFlags in) const -template +template void -MappingFE::compute_data (const UpdateFlags update_flags, - const Quadrature &q, - const unsigned int n_original_q_points, - InternalData &data) const +MappingFEField::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(); @@ -270,16 +259,15 @@ MappingFE::compute_data (const UpdateFlags update_f data.shape_second_derivatives.resize(data.n_shape_functions * n_q_points); compute_shapes_virtual (q.get_points(), data); - } -template +template void -MappingFE::compute_face_data (const UpdateFlags update_flags, - const Quadrature &q, - const unsigned int n_original_q_points, - InternalData &data) const +MappingFEField::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); @@ -341,10 +329,10 @@ MappingFE::compute_face_data (const UpdateFlags update_f } -template +template typename Mapping::InternalDataBase * -MappingFE::get_data (const UpdateFlags update_flags, - const Quadrature &quadrature) const +MappingFEField::get_data (const UpdateFlags update_flags, + const Quadrature &quadrature) const { InternalData *data = new InternalData(*fe, fe_mask); this->compute_data (update_flags, quadrature, @@ -354,10 +342,10 @@ MappingFE::get_data (const UpdateFlags update_flags, -template +template typename Mapping::InternalDataBase * -MappingFE::get_face_data (const UpdateFlags update_flags, - const Quadrature& quadrature) const +MappingFEField::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)); @@ -368,10 +356,10 @@ MappingFE::get_face_data (const UpdateFlags update_flags } -template +template typename Mapping::InternalDataBase * -MappingFE::get_subface_data (const UpdateFlags update_flags, - const Quadrature& quadrature) const +MappingFEField::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)); @@ -382,11 +370,11 @@ MappingFE::get_subface_data (const UpdateFlags update_fl } -// Note that the CellSimilarity flag is modifyable, since MappingFE can need to +// Note that the CellSimilarity flag is modifyable, since MappingFEField can need to // recalculate data even when cells are similar. -template +template void -MappingFE::fill_fe_values ( +MappingFEField::fill_fe_values ( const typename Triangulation::cell_iterator &cell, const Quadrature &q, typename Mapping::InternalDataBase &mapping_data, @@ -406,14 +394,6 @@ MappingFE::fill_fe_values ( 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; @@ -567,9 +547,9 @@ MappingFE::fill_fe_values ( -template +template void -MappingFE::fill_fe_face_values ( +MappingFEField::fill_fe_face_values ( const typename Triangulation::cell_iterator &cell, const unsigned int face_no, const Quadrature &q, @@ -580,16 +560,9 @@ MappingFE::fill_fe_face_values ( 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); @@ -608,30 +581,23 @@ MappingFE::fill_fe_face_values ( quadrature_points, JxW_values, exterior_forms, normal_vectors, jacobians, inverse_jacobians); - // quadrature_points, JxW_values, - // exterior_forms, normal_vectors); } -template +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 +MappingFEField::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 { - //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, @@ -653,16 +619,13 @@ MappingFE::fill_fe_subface_values (const typename Triang quadrature_points, JxW_values, exterior_forms, normal_vectors, jacobians, inverse_jacobians); - // quadrature_points, JxW_values, - // exterior_forms, normal_vectors); } - -template +template void -MappingFE::transform ( +MappingFEField::transform ( const VectorSlice > > input, VectorSlice > > output, const typename Mapping::InternalDataBase &mapping_data, @@ -675,9 +638,9 @@ MappingFE::transform ( -template +template void -MappingFE::transform ( +MappingFEField::transform ( const VectorSlice > > input, VectorSlice > > output, const typename Mapping::InternalDataBase &mapping_data, @@ -695,8 +658,8 @@ MappingFE::transform ( } -template -void MappingFE::transform +template +void MappingFEField::transform (const VectorSlice > > input, VectorSlice > > output, const typename Mapping::InternalDataBase &mapping_data, @@ -711,9 +674,9 @@ void MappingFE::transform } -template +template template < int rank > -void MappingFE::transform_fields( +void MappingFEField::transform_fields( const VectorSlice > > input, VectorSlice > > output, const typename Mapping::InternalDataBase &mapping_data, @@ -773,9 +736,9 @@ void MappingFE::transform_fields( } -template +template template < int rank > -void MappingFE::transform_differential_forms( +void MappingFEField::transform_differential_forms( const VectorSlice > > input, VectorSlice > > output, const typename Mapping::InternalDataBase &mapping_data, @@ -807,9 +770,9 @@ void MappingFE::transform_differential_forms( -template +template Point -MappingFE:: +MappingFEField:: transform_unit_to_real_cell (const typename Triangulation::cell_iterator &cell, const Point &p) const { @@ -819,6 +782,7 @@ transform_unit_to_real_cell (const typename Triangulation::cell_it Threads::Mutex::ScopedLock lock(mutex); update_internal_dofs(cell); + const Quadrature point_quadrature(p); std::auto_ptr mdata (dynamic_cast ( @@ -828,9 +792,9 @@ transform_unit_to_real_cell (const typename Triangulation::cell_it } -template +template Point -MappingFE:: +MappingFEField:: transform_unit_to_real_cell_internal (const InternalData &data) const { Point p_real; @@ -847,32 +811,37 @@ transform_unit_to_real_cell_internal (const InternalData &data) const -template +template Point -MappingFE:: +MappingFEField:: transform_real_to_unit_cell (const typename Triangulation::cell_iterator &cell, const Point &p) const { - { - Threads::Mutex::ScopedLock lock(mutex); - update_internal_dofs(cell); - } + 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; + try + { + initial_p_unit + = StaticMappingQ1::mapping.transform_real_to_unit_cell(cell, p); + } + catch (const typename Mapping::ExcTransformationFailed &) + { + // mirror the conditions of the code below to determine if we need to + // use an arbitrary starting point or if we just need to rethrow the + // exception + for (unsigned int d=0; d::project_to_unit_cell(initial_p_unit); + // for (unsigned int d=0; d point_quadrature(initial_p_unit); UpdateFlags update_flags = update_transformation_values|update_transformation_gradients; @@ -887,21 +856,15 @@ transform_real_to_unit_cell (const typename Triangulation::cell_it } -template +template Point -MappingFE:: +MappingFEField:: 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(); (void)n_shapes; Assert(n_shapes!=0, ExcInternalError()); @@ -911,60 +874,42 @@ transform_real_to_unit_cell_internal // 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 @@ -981,64 +926,47 @@ transform_real_to_unit_cell_internal 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; - } + 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; - } + 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 +template void -MappingFE::compute_fill ( +MappingFEField::compute_fill ( const typename Triangulation::cell_iterator &cell, const unsigned int n_q_points, const typename QProjector::DataSetDescriptor data_set, @@ -1047,8 +975,10 @@ MappingFE::compute_fill ( std::vector > &quadrature_points) const { const UpdateFlags update_flags(data.current_update_flags()); - Threads::Mutex::ScopedLock lock(mutex); - update_internal_dofs(cell); + { + Threads::Mutex::ScopedLock lock(mutex); + update_internal_dofs(cell); + } // first compute quadrature points if (update_flags & update_quadrature_points) @@ -1132,9 +1062,9 @@ MappingFE::compute_fill ( } -template +template void -MappingFE::compute_fill_face ( +MappingFEField::compute_fill_face ( const typename Triangulation::cell_iterator &cell, const unsigned int face_no, const unsigned int subface_no, @@ -1271,40 +1201,36 @@ MappingFE::compute_fill_face ( } -template +template unsigned int -MappingFE::get_degree() const +MappingFEField::get_degree() const { return fe->degree; } -template +template ComponentMask -MappingFE::get_fe_mask() const +MappingFEField::get_component_mask() const { return this->fe_mask; } -template +template Mapping * -MappingFE::clone () const +MappingFEField::clone () const { - return new MappingFE(*this); + return new MappingFEField(*this); } -template +template void -MappingFE::update_internal_dofs ( +MappingFEField::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; - } + Assert(euler_dof_handler != 0, ExcMessage("euler_dof_handler is empty")); typename DH::cell_iterator dof_cell(*cell, euler_dof_handler); Assert (dof_cell->active() == true, ExcInactiveCell()); @@ -1312,27 +1238,25 @@ MappingFE::update_internal_dofs ( dof_cell->get_dof_indices(dof_indices); for (unsigned int i=0; i -void MappingFE::update_euler_vector_using_triangulation -(VECTOR &vector) +template +void MappingFEField::update_euler_vector_using_triangulation +(VECTOR &vector) const { + AssertDimension(vector.size(), euler_dof_handler->n_dofs()); if ( fe->has_support_points() ) { - std::vector > support_points = fe->get_unit_support_points(); typename DH::active_cell_iterator cell; - Quadrature quad(support_points); + const Quadrature quad(fe->get_unit_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()); + AssertDimension(fe->dofs_per_cell, fe->get_unit_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) @@ -1343,36 +1267,45 @@ void MappingFE::update_euler_vector_using_triangulation 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]; + if (fe_mask[comp]) + vector(dofs[q]) = points[q][fe_to_real[comp]]; } } - } else { - // Construct a MappingFE with an FEQ + // Construct a MappingFEField with an FEQ, and construct a + // standard iso-parametric interpolation on the Triangulation. + // + // Once we have this, interpolate with the given finite element + // to get a Mapping which is interpolatory at the support points + // of FE_Q(fe->degree()) 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); + MappingFEField 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(); + const std::vector > &points = feq.get_unit_support_points(); - // Here construct the matrix!!!! + // Here construct the interpolation matrix from FE_Q^spacedim to + // the FiniteElement used by euler_dof_handler. + // + // The interpolation matrix is then passed to the + // VectorTools::interpolate() function to generate 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]); - } + unsigned int comp_i = fe->system_to_component_index(i).first; + if (fe_mask[comp_i]) + for (unsigned int j=0; jshape_value(i, points[j]); + } } VectorTools::interpolate(dhq, *euler_dof_handler, transfer, eulerq, vector); } @@ -1382,7 +1315,7 @@ void MappingFE::update_euler_vector_using_triangulation // explicit instantiations -#include "mapping_fe.inst" +#include "mapping_fe_field.inst" DEAL_II_NAMESPACE_CLOSE diff --git a/source/fe/mapping_fe.inst.in b/source/fe/mapping_fe_field.inst.in similarity index 82% rename from source/fe/mapping_fe.inst.in rename to source/fe/mapping_fe_field.inst.in index fbb8b6d283..ce653162ef 100644 --- a/source/fe/mapping_fe.inst.in +++ b/source/fe/mapping_fe_field.inst.in @@ -19,7 +19,7 @@ for (deal_II_dimension : DIMENSIONS; deal_II_space_dimension : SPACE_DIMENSIONS) { #if deal_II_dimension <= deal_II_space_dimension - template class MappingFE, dealii::Vector >; + template class MappingFEField, dealii::DoFHandler >; #endif } diff --git a/tests/fe/mapping_fe_field_real_to_unit_b1.cc b/tests/fe/mapping_fe_field_real_to_unit_b1.cc new file mode 100644 index 0000000000..0463f9b182 --- /dev/null +++ b/tests/fe/mapping_fe_field_real_to_unit_b1.cc @@ -0,0 +1,141 @@ +// --------------------------------------------------------------------- +// +// Copyright (C) 2006 - 2015 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 +// MappingFEField::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 +#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 feb(1); + const FESystem fesystem(feb, spacedim); + DoFHandler dhb(triangulation); + dhb.distribute_dofs(fesystem); + Vector eulerq(dhb.n_dofs()); + const ComponentMask mask(spacedim, true); + + MappingFEField map(eulerq, dhb, 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_field_real_to_unit_b1.output similarity index 100% rename from tests/fe/mapping_fe_real_to_unit_q1.output rename to tests/fe/mapping_fe_field_real_to_unit_b1.output diff --git a/tests/fe/mapping_fe_field_real_to_unit_b2_curved.cc b/tests/fe/mapping_fe_field_real_to_unit_b2_curved.cc new file mode 100644 index 0000000000..7efb2357d9 --- /dev/null +++ b/tests/fe/mapping_fe_field_real_to_unit_b2_curved.cc @@ -0,0 +1,145 @@ +// --------------------------------------------------------------------- +// +// Copyright (C) 2006 - 2015 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 +// MappingFEField::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 + +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 feb(2); + const FESystem fesystem(feb, spacedim); + DoFHandler dhb(triangulation); + dhb.distribute_dofs(fesystem); + Vector eulerq(dhb.n_dofs()); + const ComponentMask mask(spacedim, true); + + MappingFEField map(eulerq, dhb, 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_field_real_to_unit_b2_curved.output similarity index 100% rename from tests/fe/mapping_fe_real_to_unit_q5_curved.output rename to tests/fe/mapping_fe_field_real_to_unit_b2_curved.output diff --git a/tests/fe/mapping_fe_field_real_to_unit_b2_mask.cc b/tests/fe/mapping_fe_field_real_to_unit_b2_mask.cc new file mode 100644 index 0000000000..f8e4be9d7c --- /dev/null +++ b/tests/fe/mapping_fe_field_real_to_unit_b2_mask.cc @@ -0,0 +1,143 @@ +// --------------------------------------------------------------------- +// +// Copyright (C) 2006 - 2015 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 +// MappingFEField::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 +#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 fesystem(FE_Q(1),1, FE_Bernstein(2),spacedim); + + DoFHandler dhb(triangulation); + dhb.distribute_dofs(fesystem); + Vector eulerq(dhb.n_dofs()); + + // Let's use only the last spacedim components + ComponentMask mask(spacedim+1, true); + mask.set(0, false); + + MappingFEField map(eulerq, dhb, 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); + + 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<2,3>(); + + // test_real_to_unit_cell<1,3>(); + return 0; +} diff --git a/tests/fe/mapping_fe_field_real_to_unit_b2_mask.output b/tests/fe/mapping_fe_field_real_to_unit_b2_mask.output new file mode 100644 index 0000000000..5693e8cff9 --- /dev/null +++ b/tests/fe/mapping_fe_field_real_to_unit_b2_mask.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 diff --git a/tests/fe/mapping_fe_field_real_to_unit_b3_curved.cc b/tests/fe/mapping_fe_field_real_to_unit_b3_curved.cc new file mode 100644 index 0000000000..1ed55433fe --- /dev/null +++ b/tests/fe/mapping_fe_field_real_to_unit_b3_curved.cc @@ -0,0 +1,146 @@ +// --------------------------------------------------------------------- +// +// Copyright (C) 2006 - 2015 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 +// MappingFEField::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 + +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 feb(3); + const FESystem fesystem(feb, spacedim); + DoFHandler dhb(triangulation); + dhb.distribute_dofs(fesystem); + Vector eulerq(dhb.n_dofs()); + const ComponentMask mask(spacedim, true); + + MappingFEField map(eulerq, dhb, 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_field_real_to_unit_b3_curved.output b/tests/fe/mapping_fe_field_real_to_unit_b3_curved.output new file mode 100644 index 0000000000..5693e8cff9 --- /dev/null +++ b/tests/fe/mapping_fe_field_real_to_unit_b3_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 diff --git a/tests/fe/mapping_fe_field_real_to_unit_b4_curved.cc b/tests/fe/mapping_fe_field_real_to_unit_b4_curved.cc new file mode 100644 index 0000000000..2cdf1cebae --- /dev/null +++ b/tests/fe/mapping_fe_field_real_to_unit_b4_curved.cc @@ -0,0 +1,146 @@ +// --------------------------------------------------------------------- +// +// Copyright (C) 2006 - 2015 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 +// MappingFEField::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 + +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 feb(4); + const FESystem fesystem(feb, spacedim); + DoFHandler dhb(triangulation); + dhb.distribute_dofs(fesystem); + Vector eulerq(dhb.n_dofs()); + const ComponentMask mask(spacedim, true); + + MappingFEField map(eulerq, dhb, 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_field_real_to_unit_b4_curved.output b/tests/fe/mapping_fe_field_real_to_unit_b4_curved.output new file mode 100644 index 0000000000..5693e8cff9 --- /dev/null +++ b/tests/fe/mapping_fe_field_real_to_unit_b4_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 diff --git a/tests/fe/mapping_fe_field_real_to_unit_b5_curved.cc b/tests/fe/mapping_fe_field_real_to_unit_b5_curved.cc new file mode 100644 index 0000000000..7886f74919 --- /dev/null +++ b/tests/fe/mapping_fe_field_real_to_unit_b5_curved.cc @@ -0,0 +1,146 @@ +// --------------------------------------------------------------------- +// +// Copyright (C) 2006 - 2015 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 +// MappingFEField::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 + +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 feb(5); + const FESystem fesystem(feb, spacedim); + DoFHandler dhb(triangulation); + dhb.distribute_dofs(fesystem); + Vector eulerq(dhb.n_dofs()); + const ComponentMask mask(spacedim, true); + + MappingFEField map(eulerq, dhb, 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_field_real_to_unit_b5_curved.output b/tests/fe/mapping_fe_field_real_to_unit_b5_curved.output new file mode 100644 index 0000000000..5693e8cff9 --- /dev/null +++ b/tests/fe/mapping_fe_field_real_to_unit_b5_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 diff --git a/tests/fe/mapping_fe_real_to_unit_q1.cc b/tests/fe/mapping_fe_field_real_to_unit_q1.cc similarity index 95% rename from tests/fe/mapping_fe_real_to_unit_q1.cc rename to tests/fe/mapping_fe_field_real_to_unit_q1.cc index 8e0e1744a2..9ac09034dc 100644 --- a/tests/fe/mapping_fe_real_to_unit_q1.cc +++ b/tests/fe/mapping_fe_field_real_to_unit_q1.cc @@ -17,7 +17,7 @@ // 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 +// MappingFEField::transform_unit_to_real_cell that we get the same point as // we had in the beginning. #include "../tests.h" @@ -25,7 +25,7 @@ #include #include #include -#include +#include #include #include #include @@ -80,7 +80,7 @@ void test_real_to_unit_cell() Vector eulerq(dhq.n_dofs()); const ComponentMask mask(spacedim, true); - MappingFE map(eulerq, dhq, mask); + MappingFEField map(eulerq, dhq, mask); map.update_euler_vector_using_triangulation(eulerq); diff --git a/tests/fe/mapping_fe_field_real_to_unit_q1.output b/tests/fe/mapping_fe_field_real_to_unit_q1.output new file mode 100644 index 0000000000..c68c1fb460 --- /dev/null +++ b/tests/fe/mapping_fe_field_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_field_real_to_unit_q2_curved.cc b/tests/fe/mapping_fe_field_real_to_unit_q2_curved.cc new file mode 100644 index 0000000000..0cce6f055b --- /dev/null +++ b/tests/fe/mapping_fe_field_real_to_unit_q2_curved.cc @@ -0,0 +1,147 @@ +// --------------------------------------------------------------------- +// +// Copyright (C) 2006 - 2015 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 +// MappingFEField::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(2); + const FESystem fesystem(feq, spacedim); + DoFHandler dhq(triangulation); + dhq.distribute_dofs(fesystem); + Vector eulerq(dhq.n_dofs()); + const ComponentMask mask(spacedim, true); + + MappingFEField 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_field_real_to_unit_q2_curved.output b/tests/fe/mapping_fe_field_real_to_unit_q2_curved.output new file mode 100644 index 0000000000..5693e8cff9 --- /dev/null +++ b/tests/fe/mapping_fe_field_real_to_unit_q2_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 diff --git a/tests/fe/mapping_fe_field_real_to_unit_q2_mask.cc b/tests/fe/mapping_fe_field_real_to_unit_q2_mask.cc new file mode 100644 index 0000000000..bd8ef97191 --- /dev/null +++ b/tests/fe/mapping_fe_field_real_to_unit_q2_mask.cc @@ -0,0 +1,143 @@ +// --------------------------------------------------------------------- +// +// Copyright (C) 2006 - 2015 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 +// MappingFEField::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 fesystem(FE_Q(1),1, FE_Q(2),spacedim); + + DoFHandler dhq(triangulation); + dhq.distribute_dofs(fesystem); + Vector eulerq(dhq.n_dofs()); + + // Let's use only the last spacedim components + ComponentMask mask(spacedim+1, true); + mask.set(0, false); + + MappingFEField 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_field_real_to_unit_q2_mask.output b/tests/fe/mapping_fe_field_real_to_unit_q2_mask.output new file mode 100644 index 0000000000..c68c1fb460 --- /dev/null +++ b/tests/fe/mapping_fe_field_real_to_unit_q2_mask.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_field_real_to_unit_q3_curved.cc b/tests/fe/mapping_fe_field_real_to_unit_q3_curved.cc new file mode 100644 index 0000000000..e70517993f --- /dev/null +++ b/tests/fe/mapping_fe_field_real_to_unit_q3_curved.cc @@ -0,0 +1,148 @@ +// --------------------------------------------------------------------- +// +// Copyright (C) 2006 - 2015 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 +// MappingFEField::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(3); + const FESystem fesystem(feq, spacedim); + DoFHandler dhq(triangulation); + dhq.distribute_dofs(fesystem); + Vector eulerq(dhq.n_dofs()); + const ComponentMask mask(spacedim, true); + + MappingFEField 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_field_real_to_unit_q3_curved.output b/tests/fe/mapping_fe_field_real_to_unit_q3_curved.output new file mode 100644 index 0000000000..5693e8cff9 --- /dev/null +++ b/tests/fe/mapping_fe_field_real_to_unit_q3_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 diff --git a/tests/fe/mapping_fe_field_real_to_unit_q4_curved.cc b/tests/fe/mapping_fe_field_real_to_unit_q4_curved.cc new file mode 100644 index 0000000000..3ff565dfe3 --- /dev/null +++ b/tests/fe/mapping_fe_field_real_to_unit_q4_curved.cc @@ -0,0 +1,145 @@ +// --------------------------------------------------------------------- +// +// Copyright (C) 2006 - 2015 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 +// MappingFEField::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 + +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(4); + const FESystem fesystem(feq, spacedim); + DoFHandler dhq(triangulation); + dhq.distribute_dofs(fesystem); + Vector eulerq(dhq.n_dofs()); + const ComponentMask mask(spacedim, true); + + MappingFEField 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_field_real_to_unit_q4_curved.output b/tests/fe/mapping_fe_field_real_to_unit_q4_curved.output new file mode 100644 index 0000000000..5693e8cff9 --- /dev/null +++ b/tests/fe/mapping_fe_field_real_to_unit_q4_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 diff --git a/tests/fe/mapping_fe_real_to_unit_q5_curved.cc b/tests/fe/mapping_fe_field_real_to_unit_q5_curved.cc similarity index 95% rename from tests/fe/mapping_fe_real_to_unit_q5_curved.cc rename to tests/fe/mapping_fe_field_real_to_unit_q5_curved.cc index 1de8780b31..1e04f9527b 100644 --- a/tests/fe/mapping_fe_real_to_unit_q5_curved.cc +++ b/tests/fe/mapping_fe_field_real_to_unit_q5_curved.cc @@ -17,7 +17,7 @@ // 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 +// MappingFEField::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 @@ -29,7 +29,7 @@ #include #include #include -#include +#include #include #include #include @@ -37,9 +37,6 @@ #include #include -#include -#include - using namespace dealii; @@ -101,7 +98,7 @@ void test_real_to_unit_cell() Vector eulerq(dhq.n_dofs()); const ComponentMask mask(spacedim, true); - MappingFE map(eulerq, dhq, mask); + MappingFEField map(eulerq, dhq, mask); map.update_euler_vector_using_triangulation(eulerq); @@ -118,7 +115,6 @@ void test_real_to_unit_cell() // the same point again const Point 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()); diff --git a/tests/fe/mapping_fe_field_real_to_unit_q5_curved.output b/tests/fe/mapping_fe_field_real_to_unit_q5_curved.output new file mode 100644 index 0000000000..5693e8cff9 --- /dev/null +++ b/tests/fe/mapping_fe_field_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