From a306bb03411bf9759d608d69556caf2de80a4cd5 Mon Sep 17 00:00:00 2001 From: Jaeryun Yim Date: Thu, 14 Jul 2016 21:04:15 +0900 Subject: [PATCH] Add get_linear_shape(). + Edit private/protected. --- include/deal.II/fe/fe_p1nc.h | 76 ++++--- source/fe/fe_p1nc.cc | 424 +++++++++++++++-------------------- 2 files changed, 220 insertions(+), 280 deletions(-) diff --git a/include/deal.II/fe/fe_p1nc.h b/include/deal.II/fe/fe_p1nc.h index ff32ebcff6..3a95887c77 100644 --- a/include/deal.II/fe/fe_p1nc.h +++ b/include/deal.II/fe/fe_p1nc.h @@ -31,8 +31,36 @@ DEAL_II_NAMESPACE_OPEN class FE_P1NC : public FiniteElement<2,2> { -private: +public: + /** + * Constructor for nonparametric version of P1 nonconforming element. + */ + FE_P1NC() ; + + /** + * Return the name of the class for the element. + */ + virtual std::string get_name () const ; + + /** + * Return the update flags which are needed. + */ + virtual UpdateFlags requires_update_flags (const UpdateFlags flags) const ; + + /** + * Copy constructor. + */ + virtual FiniteElement<2,2> *clone () const ; + + /** + * Destructor. + */ + virtual ~FE_P1NC (); + + + +private: static std::vector @@ -47,6 +75,19 @@ private: get_dpo_vector (); + /** + * Compute the linear shape functions phi(x,y) = ax + by + c + * such that each midpoint value on two connecting edges is a half, + * and two other midpoint values are all zero. + */ + static + void + get_linear_shape (const Triangulation<2,2>::cell_iterator &cell, + std::vector &a, + std::vector &b, + std::vector &c); + + /** * Do the work which is needed before cellwise data computation. @@ -65,13 +106,13 @@ private: /** * Compute the data on the current cell. */ - virtual + virtual void fill_fe_values (const Triangulation<2,2>::cell_iterator &cell, const CellSimilarity::Similarity , const Quadrature<2> &quadrature, const Mapping<2,2> &mapping, - const Mapping<2,2>::InternalDataBase &, + const Mapping<2,2>::InternalDataBase &, const internal::FEValues::MappingRelatedData<2,2> &mapping_data, const FiniteElement<2,2>::InternalDataBase &fe_internal, internal::FEValues::FiniteElementRelatedData<2,2> &output_data) const; @@ -112,35 +153,6 @@ private: -public: - /** - * Constructor for nonparametric version of P1 nonconforming element. - */ - FE_P1NC() ; - - /** - * Return the name of the class for the element. - */ - virtual std::string get_name () const ; - - /** - * Return the update flags which are needed. - */ - virtual UpdateFlags requires_update_flags (const UpdateFlags flags) const ; - - /** - * Copy constructor. - */ - virtual FiniteElement<2,2> *clone () const ; - - /** - * Destructor. - */ - virtual ~FE_P1NC (); - - -protected: - /** * Create the constraints matrix for hanging edges. */ diff --git a/source/fe/fe_p1nc.cc b/source/fe/fe_p1nc.cc index 161b461344..e618997cdf 100644 --- a/source/fe/fe_p1nc.cc +++ b/source/fe/fe_p1nc.cc @@ -19,6 +19,61 @@ DEAL_II_NAMESPACE_OPEN +FE_P1NC::FE_P1NC() + : + FiniteElement<2,2>(FiniteElementData<2>(get_dpo_vector(), + 1, + 1, + FiniteElementData<2>::L2), + std::vector (1, false), + get_nonzero_component()) +{ + + // face support points: 2 end vertices + unit_face_support_points.resize(2); + unit_face_support_points[0][0] = 0.0 ; + unit_face_support_points[1][0] = 1.0 ; + + + // initialize constraints matrix + initialize_constraints () ; +} + + + +std::string FE_P1NC::get_name () const +{ + return "FE_P1NC"; +} + + +UpdateFlags FE_P1NC::requires_update_flags (const UpdateFlags flags) const +{ + UpdateFlags out = update_default; + + if (flags & update_values) + out |= update_values | update_quadrature_points ; + if (flags & update_gradients) + out |= update_gradients ; + if (flags & update_cell_normal_vectors) + out |= update_cell_normal_vectors | update_JxW_values; + if (flags & update_hessians) + out |= update_hessians; + if (flags & update_hessians) + out |= update_hessians; + + return out; +} + + +FiniteElement<2,2> *FE_P1NC::clone () const +{ + return new FE_P1NC(*this); +} + +FE_P1NC::~FE_P1NC () {} + + std::vector FE_P1NC::get_nonzero_component() { @@ -41,44 +96,12 @@ FE_P1NC::get_dpo_vector () } -FiniteElement<2,2>::InternalDataBase * -FE_P1NC::get_data (const UpdateFlags update_flags, - const Mapping<2,2> &, - const Quadrature<2> &, - dealii::internal::FEValues::FiniteElementRelatedData<2,2> &) const -{ - FiniteElement<2,2>::InternalDataBase *data = new FiniteElement<2,2>::InternalDataBase; - - data->update_each = requires_update_flags(update_flags); - - return data; - - - - -} - - - void -FE_P1NC::fill_fe_values (const Triangulation<2,2>::cell_iterator &cell, - const CellSimilarity::Similarity , - const Quadrature<2> &quadrature, - const Mapping<2,2> &mapping, - const Mapping<2,2>::InternalDataBase &, - const internal::FEValues::MappingRelatedData<2,2> &mapping_data, - const FiniteElement<2,2>::InternalDataBase &fe_internal, - internal::FEValues::FiniteElementRelatedData<2,2> &output_data) const +FE_P1NC::get_linear_shape (const Triangulation<2,2>::cell_iterator &cell, + std::vector &a, + std::vector &b, + std::vector &c) { - const UpdateFlags flags(fe_internal.update_each) ; - - const unsigned int n_q_points = mapping_data.quadrature_points.size(); - - std::vector values(flags & update_values ? this->dofs_per_cell : 0); - std::vector > grads(flags & update_gradients ? this->dofs_per_cell : 0); - - - // edge midpoints std::vector > mpt(4) ; @@ -100,10 +123,8 @@ FE_P1NC::fill_fe_values (const Triangulation<2,2>::cell_iterator &cell cpt(1) = (mpt[0](1) + mpt[1](1) + mpt[2](1) + mpt[3](1))/4.0 ; - // basis functions with a half value: phi(x,y) = ax + by + c - std::vector a(4), b(4), c(4) ; - double det ; + double det ; det = (mpt[0](0)-mpt[1](0))*(mpt[2](1)-mpt[3](1)) - (mpt[2](0)-mpt[3](0))*(mpt[0](1)-mpt[1](1)) ; a[0] = ((mpt[2](1)-mpt[3](1))*(0.5) -(mpt[0](1)-mpt[1](1))*(0.5))/det ; @@ -121,6 +142,51 @@ FE_P1NC::fill_fe_values (const Triangulation<2,2>::cell_iterator &cell c[2] = 0.25 - cpt(0)*a[2] - cpt(1)*b[2] ; c[3] = 0.25 - cpt(0)*a[3] - cpt(1)*b[3] ; +} + + + + +FiniteElement<2,2>::InternalDataBase * +FE_P1NC::get_data (const UpdateFlags update_flags, + const Mapping<2,2> &, + const Quadrature<2> &, + dealii::internal::FEValues::FiniteElementRelatedData<2,2> &) const +{ + FiniteElement<2,2>::InternalDataBase *data = new FiniteElement<2,2>::InternalDataBase; + + data->update_each = requires_update_flags(update_flags); + + return data; + + + + +} + + + +void +FE_P1NC::fill_fe_values (const Triangulation<2,2>::cell_iterator &cell, + const CellSimilarity::Similarity , + const Quadrature<2> &quadrature, + const Mapping<2,2> &mapping, + const Mapping<2,2>::InternalDataBase &, + const internal::FEValues::MappingRelatedData<2,2> &mapping_data, + const FiniteElement<2,2>::InternalDataBase &fe_internal, + internal::FEValues::FiniteElementRelatedData<2,2> &output_data) const +{ + const UpdateFlags flags(fe_internal.update_each) ; + + const unsigned int n_q_points = mapping_data.quadrature_points.size(); + + std::vector values(flags & update_values ? this->dofs_per_cell : 0); + std::vector > grads(flags & update_gradients ? this->dofs_per_cell : 0); + + + // linear shape with a half value: phi(x,y) = ax + by + c + std::vector a(4), b(4), c(4); + get_linear_shape (cell, a, b, c); // compute basis functions @@ -145,22 +211,22 @@ FE_P1NC::fill_fe_values (const Triangulation<2,2>::cell_iterator &cell } - // hessian - std::vector > hessians(flags & update_hessians ? this->dofs_per_cell : 0); - if (flags & update_hessians) - { - for (unsigned int i=0; i > hessians(flags & update_hessians ? this->dofs_per_cell : 0); + if (flags & update_hessians) + { + for (unsigned int i=0; idofs_per_cell; ++k) + for (unsigned int k=0; kdofs_per_cell; ++k) { - hessians[k][0][0] = 0.0 ; - hessians[k][0][1] = 0.0 ; - hessians[k][1][0] = 0.0 ; - hessians[k][1][1] = 0.0 ; - output_data.shape_hessians[k][i] = hessians[k]; + hessians[k][0][0] = 0.0 ; + hessians[k][0][1] = 0.0 ; + hessians[k][1][0] = 0.0 ; + hessians[k][1][1] = 0.0 ; + output_data.shape_hessians[k][i] = hessians[k]; } } - } + } // When this function is called for the graphic out, // MappingRelatedData does not work properly. @@ -185,13 +251,13 @@ FE_P1NC::fill_fe_values (const Triangulation<2,2>::cell_iterator &cell void FE_P1NC::fill_fe_face_values (const Triangulation<2,2>::cell_iterator &cell, - const unsigned int face_no, - const Quadrature<1> &quadrature, - const Mapping<2,2> &mapping, - const Mapping<2,2>::InternalDataBase &mapping_internal, - const dealii::internal::FEValues::MappingRelatedData<2,2> &mapping_data, - const InternalDataBase &fe_internal, - dealii::internal::FEValues::FiniteElementRelatedData<2,2> &output_data) const + const unsigned int face_no, + const Quadrature<1> &quadrature, + const Mapping<2,2> &mapping, + const Mapping<2,2>::InternalDataBase &mapping_internal, + const dealii::internal::FEValues::MappingRelatedData<2,2> &mapping_data, + const InternalDataBase &fe_internal, + dealii::internal::FEValues::FiniteElementRelatedData<2,2> &output_data) const { const UpdateFlags flags(fe_internal.update_each) ; @@ -202,49 +268,9 @@ FE_P1NC::fill_fe_face_values (const Triangulation<2,2>::cell_iterator std::vector > grads(flags & update_gradients ? this->dofs_per_cell : 0); - - // edge midpoints - std::vector > mpt(4) ; - - mpt[0](0) = (cell->vertex(0)(0) + cell->vertex(2)(0))/2.0 ; - mpt[0](1) = (cell->vertex(0)(1) + cell->vertex(2)(1))/2.0 ; - - mpt[1](0) = (cell->vertex(1)(0) + cell->vertex(3)(0))/2.0 ; - mpt[1](1) = (cell->vertex(1)(1) + cell->vertex(3)(1))/2.0 ; - - mpt[2](0) = (cell->vertex(0)(0) + cell->vertex(1)(0))/2.0 ; - mpt[2](1) = (cell->vertex(0)(1) + cell->vertex(1)(1))/2.0 ; - - mpt[3](0) = (cell->vertex(2)(0) + cell->vertex(3)(0))/2.0 ; - mpt[3](1) = (cell->vertex(2)(1) + cell->vertex(3)(1))/2.0 ; - - // center point - Point<2> cpt ; - cpt(0) = (mpt[0](0) + mpt[1](0) + mpt[2](0) + mpt[3](0))/4.0 ; - cpt(1) = (mpt[0](1) + mpt[1](1) + mpt[2](1) + mpt[3](1))/4.0 ; - - - // basis functions with a half value: phi(x,y) = ax + by + c - std::vector a(4), b(4), c(4) ; - double det ; - - det = (mpt[0](0)-mpt[1](0))*(mpt[2](1)-mpt[3](1)) - (mpt[2](0)-mpt[3](0))*(mpt[0](1)-mpt[1](1)) ; - - a[0] = ((mpt[2](1)-mpt[3](1))*(0.5) -(mpt[0](1)-mpt[1](1))*(0.5))/det ; - a[1] = ((mpt[2](1)-mpt[3](1))*(-0.5) -(mpt[0](1)-mpt[1](1))*(0.5))/det ; - a[2] = ((mpt[2](1)-mpt[3](1))*(0.5) -(mpt[0](1)-mpt[1](1))*(-0.5))/det ; - a[3] = ((mpt[2](1)-mpt[3](1))*(-0.5) -(mpt[0](1)-mpt[1](1))*(-0.5))/det ; - - b[0] = (-(mpt[2](0)-mpt[3](0))*(0.5) +(mpt[0](0)-mpt[1](0))*(0.5))/det ; - b[1] = (-(mpt[2](0)-mpt[3](0))*(-0.5) +(mpt[0](0)-mpt[1](0))*(0.5))/det ; - b[2] = (-(mpt[2](0)-mpt[3](0))*(0.5) +(mpt[0](0)-mpt[1](0))*(-0.5))/det ; - b[3] = (-(mpt[2](0)-mpt[3](0))*(-0.5) +(mpt[0](0)-mpt[1](0))*(-0.5))/det ; - - c[0] = 0.25 - cpt(0)*a[0] - cpt(1)*b[0] ; - c[1] = 0.25 - cpt(0)*a[1] - cpt(1)*b[1] ; - c[2] = 0.25 - cpt(0)*a[2] - cpt(1)*b[2] ; - c[3] = 0.25 - cpt(0)*a[3] - cpt(1)*b[3] ; - + // linear shape with a half value: phi(x,y) = ax + by + c + std::vector a(4), b(4), c(4); + get_linear_shape (cell, a, b, c); // compute basis functions @@ -276,26 +302,26 @@ FE_P1NC::fill_fe_face_values (const Triangulation<2,2>::cell_iterator if (n_q_points==0) { Quadrature<2> cellquadrature = QProjector<2>::project_to_face(quadrature, face_no) ; - for(unsigned int i=0;idofs_per_cell; ++k) - { - if (flags & update_values) - { - Point<2> realquadrature ; - - realquadrature = mapping.transform_unit_to_real_cell(cell, cellquadrature.point(i)) ; - values[k] = a[k]*realquadrature(0) + b[k]*realquadrature(1) + c[k] ; - output_data.shape_values[k][i] = values[k]; - } - - if (flags & update_gradients) - { - grads[k][0] = a[k] ; - grads[k][1] = b[k] ; - output_data.shape_gradients[k][i] = grads[k]; - } - - } + for (unsigned int i=0; idofs_per_cell; ++k) + { + if (flags & update_values) + { + Point<2> realquadrature ; + + realquadrature = mapping.transform_unit_to_real_cell(cell, cellquadrature.point(i)) ; + values[k] = a[k]*realquadrature(0) + b[k]*realquadrature(1) + c[k] ; + output_data.shape_values[k][i] = values[k]; + } + + if (flags & update_gradients) + { + grads[k][0] = a[k] ; + grads[k][1] = b[k] ; + output_data.shape_gradients[k][i] = grads[k]; + } + + } } @@ -308,14 +334,14 @@ FE_P1NC::fill_fe_face_values (const Triangulation<2,2>::cell_iterator void FE_P1NC::fill_fe_subface_values (const Triangulation<2,2>::cell_iterator &cell, - const unsigned int face_no, - const unsigned int sub_no, - const Quadrature<1> &quadrature, - const Mapping<2,2> &mapping, - const Mapping<2,2>::InternalDataBase &mapping_internal, - const dealii::internal::FEValues::MappingRelatedData<2,2> &mapping_data, - const InternalDataBase &fe_internal, - dealii::internal::FEValues::FiniteElementRelatedData<2,2> &output_data) const + const unsigned int face_no, + const unsigned int sub_no, + const Quadrature<1> &quadrature, + const Mapping<2,2> &mapping, + const Mapping<2,2>::InternalDataBase &mapping_internal, + const dealii::internal::FEValues::MappingRelatedData<2,2> &mapping_data, + const InternalDataBase &fe_internal, + dealii::internal::FEValues::FiniteElementRelatedData<2,2> &output_data) const { const UpdateFlags flags(fe_internal.update_each) ; @@ -326,49 +352,9 @@ FE_P1NC::fill_fe_subface_values (const Triangulation<2,2>::cell_iterator std::vector > grads(flags & update_gradients ? this->dofs_per_cell : 0); - - // edge midpoints - std::vector > mpt(4) ; - - mpt[0](0) = (cell->vertex(0)(0) + cell->vertex(2)(0))/2.0 ; - mpt[0](1) = (cell->vertex(0)(1) + cell->vertex(2)(1))/2.0 ; - - mpt[1](0) = (cell->vertex(1)(0) + cell->vertex(3)(0))/2.0 ; - mpt[1](1) = (cell->vertex(1)(1) + cell->vertex(3)(1))/2.0 ; - - mpt[2](0) = (cell->vertex(0)(0) + cell->vertex(1)(0))/2.0 ; - mpt[2](1) = (cell->vertex(0)(1) + cell->vertex(1)(1))/2.0 ; - - mpt[3](0) = (cell->vertex(2)(0) + cell->vertex(3)(0))/2.0 ; - mpt[3](1) = (cell->vertex(2)(1) + cell->vertex(3)(1))/2.0 ; - - // center point - Point<2> cpt ; - cpt(0) = (mpt[0](0) + mpt[1](0) + mpt[2](0) + mpt[3](0))/4.0 ; - cpt(1) = (mpt[0](1) + mpt[1](1) + mpt[2](1) + mpt[3](1))/4.0 ; - - - // basis functions with a half value: phi(x,y) = ax + by + c - std::vector a(4), b(4), c(4) ; - double det ; - - det = (mpt[0](0)-mpt[1](0))*(mpt[2](1)-mpt[3](1)) - (mpt[2](0)-mpt[3](0))*(mpt[0](1)-mpt[1](1)) ; - - a[0] = ((mpt[2](1)-mpt[3](1))*(0.5) -(mpt[0](1)-mpt[1](1))*(0.5))/det ; - a[1] = ((mpt[2](1)-mpt[3](1))*(-0.5) -(mpt[0](1)-mpt[1](1))*(0.5))/det ; - a[2] = ((mpt[2](1)-mpt[3](1))*(0.5) -(mpt[0](1)-mpt[1](1))*(-0.5))/det ; - a[3] = ((mpt[2](1)-mpt[3](1))*(-0.5) -(mpt[0](1)-mpt[1](1))*(-0.5))/det ; - - b[0] = (-(mpt[2](0)-mpt[3](0))*(0.5) +(mpt[0](0)-mpt[1](0))*(0.5))/det ; - b[1] = (-(mpt[2](0)-mpt[3](0))*(-0.5) +(mpt[0](0)-mpt[1](0))*(0.5))/det ; - b[2] = (-(mpt[2](0)-mpt[3](0))*(0.5) +(mpt[0](0)-mpt[1](0))*(-0.5))/det ; - b[3] = (-(mpt[2](0)-mpt[3](0))*(-0.5) +(mpt[0](0)-mpt[1](0))*(-0.5))/det ; - - c[0] = 0.25 - cpt(0)*a[0] - cpt(1)*b[0] ; - c[1] = 0.25 - cpt(0)*a[1] - cpt(1)*b[1] ; - c[2] = 0.25 - cpt(0)*a[2] - cpt(1)*b[2] ; - c[3] = 0.25 - cpt(0)*a[3] - cpt(1)*b[3] ; - + // linear shape with a half value: phi(x,y) = ax + by + c + std::vector a(4), b(4), c(4); + get_linear_shape (cell, a, b, c); // compute basis functions @@ -400,88 +386,30 @@ FE_P1NC::fill_fe_subface_values (const Triangulation<2,2>::cell_iterator if (n_q_points==0) { Quadrature<2> cellquadrature = QProjector<2>::project_to_subface(quadrature, face_no, sub_no) ; - for(unsigned int i=0;idofs_per_cell; ++k) - { - if (flags & update_values) - { - Point<2> realquadrature ; - - realquadrature = mapping.transform_unit_to_real_cell(cell, cellquadrature.point(i)) ; - values[k] = a[k]*realquadrature(0) + b[k]*realquadrature(1) + c[k] ; - output_data.shape_values[k][i] = values[k]; - } - - if (flags & update_gradients) - { - grads[k][0] = a[k] ; - grads[k][1] = b[k] ; - output_data.shape_gradients[k][i] = grads[k]; - } - - } - - } -} - - - -FE_P1NC::FE_P1NC() - : - FiniteElement<2,2>(FiniteElementData<2>(get_dpo_vector(), - 1, - 1, - FiniteElementData<2>::L2), - std::vector (1, false), - get_nonzero_component()) -{ - - // face support points: 2 end vertices - unit_face_support_points.resize(2); - unit_face_support_points[0][0] = 0.0 ; - unit_face_support_points[1][0] = 1.0 ; - - - // initialize constraints matrix - initialize_constraints () ; -} - - - -std::string FE_P1NC::get_name () const -{ - return "FE_P1NC"; -} + for (unsigned int i=0; idofs_per_cell; ++k) + { + if (flags & update_values) + { + Point<2> realquadrature ; + realquadrature = mapping.transform_unit_to_real_cell(cell, cellquadrature.point(i)) ; + values[k] = a[k]*realquadrature(0) + b[k]*realquadrature(1) + c[k] ; + output_data.shape_values[k][i] = values[k]; + } -UpdateFlags FE_P1NC::requires_update_flags (const UpdateFlags flags) const -{ - UpdateFlags out = update_default; - - if (flags & update_values) - out |= update_values | update_quadrature_points ; - if (flags & update_gradients) - out |= update_gradients ; - if (flags & update_cell_normal_vectors) - out |= update_cell_normal_vectors | update_JxW_values; - if (flags & update_hessians) - out |= update_hessians; - if (flags & update_hessians) - out |= update_hessians; - - return out; -} + if (flags & update_gradients) + { + grads[k][0] = a[k] ; + grads[k][1] = b[k] ; + output_data.shape_gradients[k][i] = grads[k]; + } + } -FiniteElement<2,2> *FE_P1NC::clone () const -{ - return new FE_P1NC(*this); + } } -FE_P1NC::~FE_P1NC () {} - - - void FE_P1NC::initialize_constraints () { -- 2.39.5