From: oliver Date: Tue, 25 Apr 2006 19:05:53 +0000 (+0000) Subject: Added a first implementation of the ABF elements, which are a X-Git-Tag: v8.0.0~11874 X-Git-Url: https://gitweb.dealii.org/cgi-bin/gitweb.cgi?a=commitdiff_plain;h=94981cfa1700bca29a9896498d28ca96718be5c7;p=dealii.git Added a first implementation of the ABF elements, which are a subspace of H_div. git-svn-id: https://svn.dealii.org/trunk@12894 0785d39b-7218-0410-832d-ea1e28bc413d --- diff --git a/deal.II/deal.II/source/fe/fe_abf.cc b/deal.II/deal.II/source/fe/fe_abf.cc new file mode 100644 index 0000000000..80f529a48c --- /dev/null +++ b/deal.II/deal.II/source/fe/fe_abf.cc @@ -0,0 +1,600 @@ +//--------------------------------------------------------------------------- +// $Id$ +// Version: $Name$ +// +// Copyright (C) 2003, 2004, 2005 by the deal.II authors +// +// This file is subject to QPL and may not be distributed +// without copyright and license information. Please refer +// to the file deal.II/doc/license.html for the text and +// further information on this license. +// +//--------------------------------------------------------------------------- + +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include + +#ifdef HAVE_STD_STRINGSTREAM +# include +#else +# include +#endif + +#include +using namespace std; + + +template +FE_ABF::FE_ABF (const unsigned int deg) + : + FE_PolyTensor, dim> ( + deg, + FiniteElementData(get_dpo_vector(deg), + dim, deg+1, FiniteElementData::Hdiv, 1), + std::vector(PolynomialsABF::compute_n_pols(deg), true), + std::vector >(PolynomialsABF::compute_n_pols(deg), + std::vector(dim,true))), + rt_order(deg) +{ + Assert (dim >= 2, ExcImpossibleInDim(dim)); + const unsigned int n_dofs = this->dofs_per_cell; + + this->mapping_type = this->contravariant; + // First, initialize the + // generalized support points and + // quadrature weights, since they + // are required for interpolation. + initialize_support_points(deg); + // Now compute the inverse node + //matrix, generating the correct + //basis functions from the raw + //ones. + FullMatrix M(n_dofs, n_dofs); + FETools::compute_node_matrix(M, *this); + + M.print (std::cout); + + this->inverse_node_matrix.reinit(n_dofs, n_dofs); + this->inverse_node_matrix.invert(M); + // From now on, the shape functions + // will be the correct ones, not + // the raw shape functions anymore. + + + // initialize the various matrices + for (unsigned int i=0; i::children_per_cell; ++i) + { + this->prolongation[i].reinit (n_dofs, n_dofs); + this->restriction[i].reinit (n_dofs, n_dofs); + } + + FETools::compute_embedding_matrices (*this, &this->prolongation[0]); + // initialize_restriction (); + + // TODO + std::vector > + face_embeddings(1<<(dim-1), FullMatrix(this->dofs_per_face, + this->dofs_per_face)); + //FETools::compute_face_embedding_matrices(*this, &face_embeddings[0], 0, 0); + this->interface_constraints.reinit((1<<(dim-1)) * this->dofs_per_face, + this->dofs_per_face); + unsigned int target_row=0; + for (unsigned int d=0;dinterface_constraints(target_row,j) = face_embeddings[d](i,j); + ++target_row; + } +} + + + +template +std::string +FE_ABF::get_name () const +{ + // note that the + // FETools::get_fe_from_name + // function depends on the + // particular format of the string + // this function returns, so they + // have to be kept in synch + +#ifdef HAVE_STD_STRINGSTREAM + std::ostringstream namebuf; +#else + std::ostrstream namebuf; +#endif + + namebuf << "FE_ABF<" << dim << ">(" << rt_order << ")"; + +#ifndef HAVE_STD_STRINGSTREAM + namebuf << std::ends; +#endif + return namebuf.str(); +} + + + +template +FiniteElement * +FE_ABF::clone() const +{ + return new FE_ABF(rt_order); +} + + +//--------------------------------------------------------------------------- +// Auxiliary and internal functions +//--------------------------------------------------------------------------- + + +#if deal_II_dimension == 1 + +template +void +FE_ABF::initialize_support_points (const unsigned int deg) +{ + return; + + Assert (false, ExcNotImplemented()); + + QGauss cell_quadrature(deg+1); + const unsigned int n_interior_points + = (deg>0) ? cell_quadrature.n_quadrature_points : 0; + + this->generalized_support_points.resize (2 + n_interior_points); + + // Number of the point being entered + unsigned int current = 0; + + + if (deg==0) return; + + interior_weights.reinit(TableIndices<3>(2+n_interior_points, 0, dim)); + + for (unsigned int k=0;kgeneralized_support_points[current++] = cell_quadrature.point(k); + + Assert (current == this->generalized_support_points.size(), + ExcInternalError()); +} + +#else + +// Version for 2d and higher. See above for 1d version +template +void +FE_ABF::initialize_support_points (const unsigned int deg) +{ + QGauss cell_quadrature(deg+2); + const unsigned int n_interior_points = cell_quadrature.n_quadrature_points; + + unsigned int n_face_points = (dim>1) ? 1 : 0; + // compute (deg+1)^(dim-1) + for (unsigned int d=1;dgeneralized_support_points.resize (GeometryInfo::faces_per_cell*n_face_points + + n_interior_points); + this->generalized_face_support_points.resize (n_face_points); + + + // These might be required when the faces contribution is computed + // Therefore they will be initialised at this point. + std::vector* > polynomials_abf(dim); + + // Generate x_1^{i} x_2^{r+1} ... + for (unsigned int dd=0; dd > > poly(dim); + for (unsigned int d=0;d (deg+1)); + poly[dd] = Polynomials::Monomial::generate_complete_basis(deg); + + polynomials_abf[dd] = new AnisotropicPolynomials(poly); + } + + // Number of the point being entered + unsigned int current = 0; + + if (dim>1) + { + QGauss face_points (deg+1); + TensorProductPolynomials legendre + = Polynomials::Legendre::generate_complete_basis(deg); + + boundary_weights.reinit(n_face_points, legendre.n()); + +// Assert (face_points.n_quadrature_points == this->dofs_per_face, +// ExcInternalError()); + + for (unsigned int k=0;kgeneralized_face_support_points[k] = face_points.point(k); + // Compute its quadrature + // contribution for each + // moment. + for (unsigned int i=0;i faces = QProjector::project_to_all_faces(face_points); + for (;current::faces_per_cell*n_face_points; + ++current) + { + // Enter the support point + // into the vector + this->generalized_support_points[current] = faces.point(current); + } + + + // Now initialise edge interior weights for the ABF elements. + // These are completely independent from the usual edge moments. They + // stem from applying the Gauss theorem to the nodal values, which + // was necessary to cast the ABF elements into the deal.II framework + // for vector valued elements. + boundary_weights_abf.reinit(faces.n_quadrature_points, polynomials_abf[0]->n() * dim); + for (unsigned int k=0;k < faces.n_quadrature_points;++k) + { + for (unsigned int i=0;in() * dim;++i) + { + boundary_weights_abf(k,i) = polynomials_abf[i%dim]-> + compute_value(i / dim, faces.point(k)) * faces.weight(k); + } + } + } + + // Create Legendre basis for the + // space D_xi Q_k + if (deg>0) + { + std::vector* > polynomials(dim); + + for (unsigned int dd=0;dd > > poly(dim); + for (unsigned int d=0;d(poly); + } + + interior_weights.reinit(TableIndices<3>(n_interior_points, polynomials[0]->n(), dim)); + + for (unsigned int k=0;kn();++i) + for (unsigned int d=0;dcompute_value(i,cell_quadrature.point(k)); + } + + for (unsigned int d=0;dgeneralized_support_points[current++] = cell_quadrature.point(k); + + // Additional functionality for the ABF elements + // TODO: Here the canonical extension of the principle + // behind the ABF elements is implemented. It is unclear, + // if this really leads to the ABF spaces in 3D! + interior_weights_abf.reinit(TableIndices<3>(cell_quadrature.n_quadrature_points, + polynomials_abf[0]->n() * dim, dim)); + Tensor<1, dim> poly_grad; + + for (unsigned int k=0;kn() * dim;++i) + { + poly_grad = polynomials_abf[i%dim]->compute_grad(i / dim,cell_quadrature.point(k)) + * cell_quadrature.weight(k); + // The minus sign comes from the use of the Gauss theorem to replace the divergence. + for (unsigned int d=0;dgeneralized_support_points.size(), + ExcInternalError()); +} + +#endif + + +#if deal_II_dimension == 1 + +template <> +std::vector +FE_ABF<1>::get_dpo_vector (const unsigned int) +{ + Assert (false, ExcImpossibleInDim(1)); + return std::vector(); +} + +#endif + + +template +std::vector +FE_ABF::get_dpo_vector (const unsigned int rt_order) +{ + // the element is face-based (not + // to be confused with George + // W. Bush's Faith Based + // Initiative...), and we have + // (rt_order+1)^(dim-1) DoFs per face + unsigned int dofs_per_face = 1; + for (unsigned int d=0; d dpo(dim+1); + dpo[dim-1] = dofs_per_face; + dpo[dim] = interior_dofs; + + return dpo; +} + + + +template +UpdateFlags +FE_ABF::update_once (const UpdateFlags) const +{ + // even the values have to be + // computed on the real cell, so + // nothing can be done in advance + return update_default; +} + + + +template +UpdateFlags +FE_ABF::update_each (const UpdateFlags flags) const +{ + UpdateFlags out = update_default; + + + if (flags & update_values) + out |= update_values | update_covariant_transformation + | update_contravariant_transformation + | update_JxW_values; + if (flags & update_gradients) + out |= update_gradients | update_covariant_transformation + | update_contravariant_transformation + | update_JxW_values; + //TODO: Set update flags appropriately and figure out, how the second + // derivatives for the RT elements can be computed correctly. + if (flags & update_second_derivatives) + out |= update_second_derivatives | update_contravariant_transformation; + + return out; +} + +//--------------------------------------------------------------------------- +// Data field initialization +//--------------------------------------------------------------------------- + + + + +template +unsigned int +FE_ABF::n_base_elements () const +{ + return 1; +} + + + +template +const FiniteElement & +FE_ABF::base_element (const unsigned int index) const +{ + Assert (index==0, ExcIndexRange(index, 0, 1)); + return *this; +} + + + +template +unsigned int +FE_ABF::element_multiplicity (const unsigned int index) const +{ + Assert (index==0, ExcIndexRange(index, 0, 1)); + return 1; +} + + + +template +bool +FE_ABF::has_support_on_face (const unsigned int shape_index, + const unsigned int face_index) const +{ + Assert (shape_index < this->dofs_per_cell, + ExcIndexRange (shape_index, 0, this->dofs_per_cell)); + Assert (face_index < GeometryInfo::faces_per_cell, + ExcIndexRange (face_index, 0, GeometryInfo::faces_per_cell)); + + // Return computed values if we + // know them easily. Otherwise, it + // is always safe to return true. + switch (rt_order) + { + case 0: + { + switch (dim) + { + case 2: + { + // only on the one + // non-adjacent face + // are the values + // actually zero. list + // these in a table + return (face_index != GeometryInfo::opposite_face[shape_index]); + } + + default: + return true; + }; + }; + + default: // other rt_order + return true; + }; + + return true; +} + + + +template +void +FE_ABF::interpolate( + std::vector&, + const std::vector&) const +{ + Assert(false, ExcNotImplemented()); +} + + + +template +void +FE_ABF::interpolate( + std::vector& local_dofs, + const std::vector >& values, + unsigned int offset) const +{ + Assert (values.size() == this->generalized_support_points.size(), + ExcDimensionMismatch(values.size(), this->generalized_support_points.size())); + Assert (local_dofs.size() == this->dofs_per_cell, + ExcDimensionMismatch(local_dofs.size(),this->dofs_per_cell)); + Assert (values[0].size() >= offset+this->n_components(), + ExcDimensionMismatch(values[0].size(),offset+this->n_components())); + + std::fill(local_dofs.begin(), local_dofs.end(), 0.); + + const unsigned int n_face_points = boundary_weights.size(0); + for (unsigned int face=0;face::faces_per_cell;++face) + for (unsigned int k=0;kdofs_per_face] += boundary_weights(k,i) + * values[face*n_face_points+k](GeometryInfo::unit_normal_direction[face]+offset); + } + + const unsigned start_cell_dofs = GeometryInfo::faces_per_cell*this->dofs_per_face; + const unsigned start_cell_points = GeometryInfo::faces_per_cell*n_face_points; + + for (unsigned int k=0;k +void +FE_ABF::interpolate( + std::vector& local_dofs, + const VectorSlice > >& values) const +{ + Assert (values.size() == this->n_components(), + ExcDimensionMismatch(values.size(), this->n_components())); + Assert (values[0].size() == this->generalized_support_points.size(), + ExcDimensionMismatch(values[0].size(), this->generalized_support_points.size())); + Assert (local_dofs.size() == this->dofs_per_cell, + ExcDimensionMismatch(local_dofs.size(),this->dofs_per_cell)); + + std::fill(local_dofs.begin(), local_dofs.end(), 0.); + + const unsigned int n_face_points = boundary_weights.size(0); + for (unsigned int face=0;face::faces_per_cell;++face) + for (unsigned int k=0;kdofs_per_face] += boundary_weights(k,i) + * values[GeometryInfo::unit_normal_direction[face]][face*n_face_points+k]; + } + + const unsigned start_cell_dofs = GeometryInfo::faces_per_cell*this->dofs_per_face; + const unsigned start_cell_points = GeometryInfo::faces_per_cell*n_face_points; + + for (unsigned int k=0;k::faces_per_cell; ++face) + { + double n_orient = (double) GeometryInfo::unit_normal_orientation[face]; + for (unsigned int fp=0; fp < n_face_points; ++fp) + { + // TODO: Check what the face_orientation has to be in 3D + unsigned int k = QProjector::DataSetDescriptor::face (face, false, n_face_points); + for (unsigned int i=0; i::unit_normal_direction[face]][k + fp]; + } + } + + // TODO: Check if this "correction" can be removed. + for (unsigned int i=0; i +unsigned int +FE_ABF::memory_consumption () const +{ + Assert (false, ExcNotImplemented ()); + return 0; +} + + +template class FE_ABF;