From: Jaeryun Yim Date: Thu, 5 Nov 2015 02:46:07 +0000 (+0900) Subject: Add a first version of the P1 non-conforming element. X-Git-Tag: v8.5.0-rc1~624^2~42 X-Git-Url: https://gitweb.dealii.org/cgi-bin/gitweb.cgi?a=commitdiff_plain;h=ed3d6f43880b4fc781ce98ef66a319624adee7f0;p=dealii.git Add a first version of the P1 non-conforming element. This version currently only covers the parameteric version of the element. --- diff --git a/include/deal.II/fe/fe_p1nc.h b/include/deal.II/fe/fe_p1nc.h new file mode 100644 index 0000000000..d4d39fe16a --- /dev/null +++ b/include/deal.II/fe/fe_p1nc.h @@ -0,0 +1,156 @@ +// --------------------------------------------------------------------- +// +// Copyright (C) 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. +// +// --------------------------------------------------------------------- + +#ifndef dealii__fe_p1nc_h +#define dealii__fe_p1nc_h + +#include +#include +#include +#include + + +DEAL_II_NAMESPACE_OPEN + + +/*!@addtogroup fe */ +/*@{*/ + + +// parametric version of P1 Nonconforming FE +class FE_P1NCParametric : public FiniteElement<2,2> +{ +private: + + + static + std::vector + get_nonzero_component(); + + + static + std::vector + get_dpo_vector (); + + + class InternalData : public FiniteElement<2,2>::InternalDataBase + { + public: + + std::vector > shape_values; + std::vector > > shape_gradients; + }; + + + virtual FiniteElement<2,2>::InternalDataBase * + get_data (const UpdateFlags update_flags, + const Mapping<2,2> &mapping, + const Quadrature<2> &quadrature) const ; + + + + + virtual + void + fill_fe_values (const Mapping<2,2> &mapping, + const Triangulation<2,2>::cell_iterator &cell, + const Quadrature<2> &quadrature, + const Mapping<2,2>::InternalDataBase &mapping_internal, + const FiniteElement<2,2>::InternalDataBase &fe_internal, + const internal::FEValues::MappingRelatedData<2,2> &mapping_data, + internal::FEValues::FiniteElementRelatedData<2,2> &output_data, + const CellSimilarity::Similarity cell_similarity) const; + + + + + virtual + void + fill_fe_face_values (const Mapping<2,2> &mapping, + const Triangulation<2,2>::cell_iterator &cell, + const unsigned int face_no, + const Quadrature<1> &quadrature, + const Mapping<2,2>::InternalDataBase &mapping_internal, + const FiniteElement<2,2>::InternalDataBase &fe_internal, + const internal::FEValues::MappingRelatedData<2,2> &mapping_data, + internal::FEValues::FiniteElementRelatedData<2,2> &output_data) const; + + + + + virtual + void + fill_fe_subface_values (const Mapping<2,2> &mapping, + 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>::InternalDataBase &mapping_internal, + const FiniteElement<2,2>::InternalDataBase &fe_internal, + const internal::FEValues::MappingRelatedData<2,2> &mapping_data, + internal::FEValues::FiniteElementRelatedData<2,2> &output_data) const; + + + +public: + FE_P1NCParametric() ; + + + + virtual std::string get_name () const ; + + + + virtual UpdateFlags update_once (const UpdateFlags flags) const ; + + virtual UpdateFlags update_each (const UpdateFlags flags) const ; + + const std::vector > &get_unit_face_support_points () const ; + + Point<1> unit_face_support_point (const unsigned index) const ; + + bool has_face_support_points () const ; + + virtual FiniteElement<2,2> *clone () const ; + + virtual ~FE_P1NCParametric (); + + + + + + virtual double shape_value (const unsigned int i, + const Point<2> &p) const ; + + virtual Tensor<1,2> shape_grad (const unsigned int i, + const Point<2> &p) const ; + + + +protected: + + void initialize_constraints () ; + +}; + + + + + +/** @}*/ + +DEAL_II_NAMESPACE_CLOSE + +#endif diff --git a/source/fe/CMakeLists.txt b/source/fe/CMakeLists.txt index 0e8f549797..e0542b863b 100644 --- a/source/fe/CMakeLists.txt +++ b/source/fe/CMakeLists.txt @@ -33,6 +33,7 @@ SET(_src fe_nothing.cc fe_poly.cc fe_poly_tensor.cc + fe_p1nc.cc fe_q_base.cc fe_q.cc fe_q_bubbles.cc diff --git a/source/fe/fe_p1nc.cc b/source/fe/fe_p1nc.cc new file mode 100644 index 0000000000..c20bdb5c63 --- /dev/null +++ b/source/fe/fe_p1nc.cc @@ -0,0 +1,477 @@ +// --------------------------------------------------------------------- +// +// Copyright (C) 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. +// +// --------------------------------------------------------------------- + + +#include + +DEAL_II_NAMESPACE_OPEN + + +// FE_P1NCParametric + +std::vector +FE_P1NCParametric::get_nonzero_component() +{ + const unsigned int dofs_per_cell = 4; + std::vector masks (dofs_per_cell); + for (unsigned int i=0; i +FE_P1NCParametric::get_dpo_vector () +{ + std::vector dpo(3); + dpo[0] = 1; //dofs per object: vertex + dpo[1] = 0; // line + dpo[2] = 0; // quad + return dpo; +} + + +FiniteElement<2,2>::InternalDataBase * +FE_P1NCParametric::get_data (const UpdateFlags update_flags, + const Mapping<2,2> &mapping, + const Quadrature<2> &quadrature) const +{ + + InternalData *data = new InternalData; + + + data->update_once = update_once(update_flags); + data->update_each = update_each(update_flags); + data->update_flags = data->update_once | data->update_each; + + const UpdateFlags flags(data->update_flags); + const unsigned int n_q_points = quadrature.size(); + + + // initialize shape value + if (flags & update_values) + data->shape_values.resize (this->dofs_per_cell, + std::vector (n_q_points)); + + // initialize shape gradient + if (flags & update_gradients) + data->shape_gradients.resize (this->dofs_per_cell, + std::vector > (n_q_points)); + + + // update + if (flags & (update_values | update_gradients)) + for (unsigned int i=0; idofs_per_cell; ++k) + { + data->shape_values[k][i] = shape_value(k,quadrature.point(i)); + } + + + if (flags & update_gradients) + for (unsigned int k=0; kdofs_per_cell; ++k) + { + data->shape_gradients[k][i] = shape_grad(k,quadrature.point(i)); + } + + } + return data; + +} + + +// get fe value in real cell from fedata in reference cell +void +FE_P1NCParametric::fill_fe_values (const Mapping<2,2> &mapping, + const Triangulation<2,2>::cell_iterator &cell, + const Quadrature<2> &quadrature, + const Mapping<2,2>::InternalDataBase &mapping_internal, + const FiniteElement<2,2>::InternalDataBase &fe_internal, + const internal::FEValues::MappingRelatedData<2,2> &mapping_data, + internal::FEValues::FiniteElementRelatedData<2,2> &output_data, + const CellSimilarity::Similarity cell_similarity) const + +{ + + Assert (dynamic_cast (&fe_internal) != 0, ExcInternalError()); + const InternalData &fe_data = static_cast (fe_internal); + + const UpdateFlags flags(fe_data.current_update_flags()); + + // gets shape function values from its input argument fedata + for (unsigned int k=0; kdofs_per_cell; ++k) + { + // shape value + if (flags & update_values) + for (unsigned int i=0; i &mapping, + const Triangulation<2,2>::cell_iterator &cell, + const unsigned int face_no, + const Quadrature<1> &quadrature, + const Mapping<2,2>::InternalDataBase &mapping_internal, + const FiniteElement<2,2>::InternalDataBase &fe_internal, + const internal::FEValues::MappingRelatedData<2,2> &mapping_data, + internal::FEValues::FiniteElementRelatedData<2,2> &output_data) const + + +{ + + Assert (dynamic_cast (&fe_internal) != 0, ExcInternalError()); + const InternalData &fe_data = static_cast (fe_internal); + + + const QProjector<2>::DataSetDescriptor offset + = QProjector<2>::DataSetDescriptor::face (face_no, + cell->face_orientation(face_no), + cell->face_flip(face_no), + cell->face_rotation(face_no), + quadrature.size()); + + const UpdateFlags flags(fe_data.update_once | fe_data.update_each); + + for (unsigned int k=0; kdofs_per_cell; ++k) + { + if (flags & update_values) + for (unsigned int i=0; i &mapping, + 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>::InternalDataBase &mapping_internal, + const FiniteElement<2,2>::InternalDataBase &fe_internal, + const internal::FEValues::MappingRelatedData<2,2> &mapping_data, + internal::FEValues::FiniteElementRelatedData<2,2> &output_data) const + +{ + + Assert (dynamic_cast (&fe_internal) != 0, ExcInternalError()); + const InternalData &fe_data = static_cast (fe_internal); + + + + const QProjector<2>::DataSetDescriptor offset + = QProjector<2>::DataSetDescriptor::subface (face_no, sub_no, + cell->face_orientation(face_no), + cell->face_flip(face_no), + cell->face_rotation(face_no), + quadrature.size(), + cell->subface_case(face_no)); + + const UpdateFlags flags(fe_data.update_once | fe_data.update_each); + + for (unsigned int k=0; kdofs_per_cell; ++k) + { + if (flags & update_values) + for (unsigned int i=0; i(FiniteElementData<2>(get_dpo_vector(), + 1, + 1, + FiniteElementData<2>::L2, + numbers::invalid_unsigned_int), + std::vector (1, false), + get_nonzero_component()) +{ + + // support points: 4 vertices + unit_support_points.resize(4) ; + unit_support_points[0][0] = 0.0 ; + unit_support_points[0][1] = 0.0 ; + unit_support_points[1][0] = 1.0 ; + unit_support_points[1][1] = 0.0 ; + unit_support_points[2][0] = 0.0 ; + unit_support_points[2][1] = 1.0 ; + unit_support_points[3][0] = 1.0 ; + unit_support_points[3][1] = 1.0 ; + + // 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 ; + + + // CONSTRAINTS MATRIX INITIALIZATION + initialize_constraints () ; +} + + + +std::string FE_P1NCParametric::get_name () const +{ + return "FE_P1NCParametric"; +} + + + +UpdateFlags FE_P1NCParametric::update_once (const UpdateFlags flags) const +{ + return (update_default | (flags & update_values)); +} + +UpdateFlags FE_P1NCParametric::update_each (const UpdateFlags flags) const +{ + UpdateFlags out = update_default; + + if (flags & update_gradients) + out |= update_gradients | update_covariant_transformation; + if (flags & update_cell_normal_vectors) + out |= update_cell_normal_vectors | update_JxW_values; + + return out; +} + + +const std::vector > &FE_P1NCParametric::get_unit_face_support_points () const +{ + + Assert ((unit_face_support_points.size() == 0) || + (unit_face_support_points.size() == this->dofs_per_face), + ExcInternalError()); + return unit_face_support_points; +} + +Point<1> FE_P1NCParametric::unit_face_support_point (const unsigned index) const +{ + Assert (index < this->dofs_per_face, + ExcIndexRange (index, 0, this->dofs_per_face)); + Assert (unit_face_support_points.size() == this->dofs_per_face, + ExcFEHasNoSupportPoints ()); + return unit_face_support_points[index]; +} + +bool FE_P1NCParametric::has_face_support_points () const +{ + return (unit_face_support_points.size() != 0); +} + +FiniteElement<2,2> *FE_P1NCParametric::clone () const +{ + return new FE_P1NCParametric(*this); +} + +FE_P1NCParametric::~FE_P1NCParametric () {} + +// Define the exact value of the shape function for the P1NC scalar +// function in the unit cell = (0,1)^2. +// +// 2-----------------|------------------3 +// | | +// | | +// | | +// | | +// | | +// | | +// | | +// | | +// | | +// - - +// | | +// | | +// | | +// | | +// | | +// | | +// | | +// | | +// | | +// 0-----------------|------------------1 +// + + +double FE_P1NCParametric::shape_value (const unsigned int i, + const Point<2> &p) const +{ + + + double value = 0.0 ; + + // P1NC with HALF COEFFICIENT + switch (i) + { + case 0: + value = -p[0]-p[1]+1.5 ; + break ; + case 1: + value = p[0]-p[1]+0.5 ; + break ; + case 2: + value = -p[0]+p[1]+0.5 ; + break ; + case 3: + value = p[0]+p[1]-0.5 ; + break ; + default: + value = 0.0 ; + }; + + return value*0.5 ; + + + // // Q1 for test + // switch (i) + // { + // case 0: + // return (1.0-p[0])*(1.0-p[1]) ; + // case 1: + // return p[0]*(1.0-p[1]) ; + // case 2: + // return (1.0-p[0])*p[1] ; + // case 3: + // return p[0]*p[1] ; + // default: + // return 0 ; + // }; + + +} + + + + +Tensor<1,2> FE_P1NCParametric::shape_grad (const unsigned int i, + const Point<2> &p) const +{ + Tensor<1,2> grad ; + + + // P1NC with HALF COEFFICIENT + switch (i) + { + case 0: + grad[0] = -1.0 ; + grad[1] = -1.0 ; + break ; + case 1: + grad[0] = 1.0 ; + grad[1] = -1.0 ; + break ; + case 2: + grad[0] = -1.0 ; + grad[1] = 1.0 ; + break ; + case 3: + grad[0] = 1.0 ; + grad[1] = 1.0 ; + break ; + default: + grad = 0.0 ; + }; + + return grad*0.5 ; + + + // // Q1 for test + // switch (i) + // { + // case 0: + // grad[0] = p[1]-1.0 ; + // grad[1] = p[0]-1.0 ; + // return grad ; + // case 1: + // grad[0] = 1.0-p[1] ; + // grad[1] = -p[0] ; + // return grad ; + // case 2: + // grad[0] = -p[1] ; + // grad[1] = 1.0-p[0] ; + // return grad ; + // case 3: + // grad[0] = p[1] ; + // grad[1] = p[0] ; + // return grad ; + // default: + // return grad ; + // }; + + +} + + +// 2015 05 13 CONSTRAINTS MATRIX FOR HANGING NODE +// BASED ON 'fe_q_base.cc' + +void FE_P1NCParametric::initialize_constraints () +{ + + std::vector > constraint_points; + // Add midpoint + constraint_points.push_back (Point<1> (0.5)); + + // Now construct relation between destination (child) and source (mother) + // dofs. + + interface_constraints + .TableBase<2,double>::reinit (interface_constraints_size()); + + + interface_constraints(0,0) = 0.5 ; + interface_constraints(0,1) = 0.5 ; + +} + + + +DEAL_II_NAMESPACE_CLOSE