From: wolf Date: Mon, 10 Jun 2002 10:39:55 +0000 (+0000) Subject: Merge the changes made to the branch named branch_non_trivial_fes up to the tag branc... X-Git-Url: https://gitweb.dealii.org/cgi-bin/gitweb.cgi?a=commitdiff_plain;h=a590830897db6d3e8aaf2da0d7ac86dcd806348e;p=dealii-svn.git Merge the changes made to the branch named branch_non_trivial_fes up to the tag branch_non_trivial_fes_mergepoints1. For the commit messages of individual changes, see the messages on this branch. git-svn-id: https://svn.dealii.org/trunk@6040 0785d39b-7218-0410-832d-ea1e28bc413d --- diff --git a/deal.II/deal.II/include/dofs/dof_accessor.templates.h b/deal.II/deal.II/include/dofs/dof_accessor.templates.h index 87c6f6ec88..2607bb7c0e 100644 --- a/deal.II/deal.II/include/dofs/dof_accessor.templates.h +++ b/deal.II/deal.II/include/dofs/dof_accessor.templates.h @@ -165,7 +165,7 @@ DoFObjectAccessor<1,dim>::get_dof_indices (std::vector &dof_indice // freedom on sub-objects which are // not allocated for this // non-active thing - Assert (this->active() || + Assert (!this->has_children() || (this->dof_handler->get_fe().dofs_per_cell == 2*this->dof_handler->get_fe().dofs_per_vertex), typename DoFAccessor::ExcNotActive()); @@ -278,7 +278,7 @@ DoFObjectAccessor<2,dim>::get_dof_indices (std::vector &dof_indice // freedom on sub-objects which are // not allocated for this // non-active thing - Assert (this->active() || + Assert (!this->has_children() || (this->dof_handler->get_fe().dofs_per_cell == 4*this->dof_handler->get_fe().dofs_per_vertex), typename DoFAccessor::ExcNotActive()); @@ -410,7 +410,7 @@ DoFObjectAccessor<3,dim>::get_dof_indices (std::vector &dof_indice // freedom on sub-objects which are // not allocated for this // non-active thing - Assert (this->active() || + Assert (!this->has_children() || (this->dof_handler->get_fe().dofs_per_cell == 8*this->dof_handler->get_fe().dofs_per_vertex), typename DoFAccessor::ExcNotActive()); diff --git a/deal.II/deal.II/include/fe/fe.h b/deal.II/deal.II/include/fe/fe.h index 6b67093942..467bb0e713 100644 --- a/deal.II/deal.II/include/fe/fe.h +++ b/deal.II/deal.II/include/fe/fe.h @@ -1,4 +1,3 @@ - //---------------------------- fe.h --------------------------- // $Id$ // Version: $Name$ @@ -65,7 +64,8 @@ class FiniteElement : public FiniteElementBase * Constructor */ FiniteElement (const FiniteElementData &fe_data, - const std::vector &restriction_is_additive_flags); + const std::vector &restriction_is_additive_flags, + const std::vector > &nonzero_components); /** * Virtual destructor. Makes sure @@ -152,16 +152,37 @@ class FiniteElement : public FiniteElementBase protected: /** - * Compute flags for initial - * update only. - * @ref{FEValuesBase} + * Given a set of flags indicating + * what quantities are requested + * from a @p{FEValues} object, + * return which of these can be + * precomputed once and for + * all. Often, the values of + * shape function at quadrature + * points can be precomputed, for + * example, in which case the + * return value of this function + * would be the logical and of + * the input @p{flags} and + * @p{update_values}. */ virtual UpdateFlags update_once (const UpdateFlags flags) const = 0; /** - * Compute flags for update on - * each cell. - * @ref{FEValuesBase} + * This is the opposite to the + * above function: given a set of + * flags indicating what we want + * to know, return which of these + * need to be computed each time + * we visit a new cell. + * + * If for the computation of one + * quantity something else is + * also required (for example, we + * often need the covariant + * transformation when gradients + * need to be computed), include + * this in the result as well. */ virtual UpdateFlags update_each (const UpdateFlags flags) const = 0; diff --git a/deal.II/deal.II/include/fe/fe_base.h b/deal.II/deal.II/include/fe/fe_base.h index cb0f732e56..7d1c32fe57 100644 --- a/deal.II/deal.II/include/fe/fe_base.h +++ b/deal.II/deal.II/include/fe/fe_base.h @@ -25,6 +25,8 @@ template class FESystem; + + /** * Dimension independent data for finite elements. See the derived * class @ref{FiniteElementBase} class for information on its use. All @@ -117,11 +119,13 @@ class FiniteElementData const unsigned int components; /** - * Default constructor. Constructs - * an element - * which is not so useful. Checking - * @p{dofs_per_cell} is therefore a good way to - * check if something went wrong. + * Default + * constructor. Constructs an + * element which is not so + * useful. Checking + * @p{dofs_per_cell} is therefore + * a good way to check if + * something went wrong. */ FiniteElementData (); @@ -179,11 +183,8 @@ class FiniteElementData unsigned int n_components () const; /** - * Comparison operator. It is not clear to - * me (WB) why we have to declare and implement - * this one explicitly. + * Comparison operator. */ -//TODO:[WB] (compiler) remove operator and let the compiler generate it as soon as it is willing to do so bool operator == (const FiniteElementData &) const; /** @@ -350,7 +351,7 @@ class FiniteElementData * introduced from the two sides are unique; it is able to handle the fact * that the constraints for some of the dofs are entered more than once. * - * @author Wolfgang Bangerth, 1998, Ralf Hartmann, Guido Kanschat, 2001 + * @author Wolfgang Bangerth, 1998, 2002, Ralf Hartmann, Guido Kanschat, 2001 */ template class FiniteElementBase : public Subscriptor, @@ -358,10 +359,15 @@ class FiniteElementBase : public Subscriptor, { public: /** - * Basis class for internal data. + * Base class for internal data. * Adds data for second derivatives to * @ref{Mapping::InternalDataBase} * + * For information about the + * general purpose of this class, + * see the documentation of the + * base class. + * * @author Guido Kanschat, 2001 */ class InternalDataBase : public Mapping::InternalDataBase @@ -402,30 +408,65 @@ class FiniteElementBase : public Subscriptor, }; /** - * Construct an object of this type. - * You have to set the - * matrices explicitly after calling - * this base class' constructor. + * Construct an object of this + * type. You have to set some + * member variables, for example + * some matrices, explicitly + * after calling this base class' + * constructor. */ FiniteElementBase (const FiniteElementData &fe_data, - const std::vector &restriction_is_additive_flags); + const std::vector &restriction_is_additive_flags, + const std::vector > &nonzero_components); /** * Return the value of the * @p{i}th shape function at the - * point @p{p}. @p{p} is a point - * on the reference element. + * point @p{p}. @p{p} is a point + * on the reference element. If + * the finite element is + * vector-valued, then return the + * value of the only non-zero + * component of the vector value + * of this shape function. If the + * shape function has more than + * one non-zero component (which + * we refer to with the term + * non-primitive), then derived + * classes implementing this + * function should throw an + * exception of type + * @p{ExcShapeFunctionNotPrimitive}. In + * that case, use the + * @ref{shape_value_component} + * function. * * An * @p{ExcUnitShapeValuesDoNotExist} * is thrown if the shape values * of the @p{FiniteElement} under - * consideration depend on the + * consideration depends on the * shape of the cell in real * space. */ - virtual double shape_value (const unsigned int i, - const Point &p) const; + virtual double shape_value (const unsigned int i, + const Point &p) const; + + /** + * Just like for @p{shape_value}, + * but this function will be + * called when the shape function + * has more than one non-zero + * vector component. In that + * case, this function should + * return the value of the + * @p{component}-th vector + * component of the @p{i}th shape + * function at point @p{p}. + */ + virtual double shape_value_component (const unsigned int i, + const Point &p, + const unsigned int component) const; /** * Return the gradient of the @@ -435,19 +476,51 @@ class FiniteElementBase : public Subscriptor, * likewise the gradient is the * gradient on the unit cell with * respect to unit cell - * coordinates. + * coordinates. If + * the finite element is + * vector-valued, then return the + * value of the only non-zero + * component of the vector value + * of this shape function. If the + * shape function has more than + * one non-zero component (which + * we refer to with the term + * non-primitive), then derived + * classes implementing this + * function should throw an + * exception of type + * @p{ExcShapeFunctionNotPrimitive}. In + * that case, use the + * @ref{shape_grad_component} + * function. * * An * @p{ExcUnitShapeValuesDoNotExist} * is thrown if the shape values * of the @p{FiniteElement} under - * consideration depend on the + * consideration depends on the * shape of the cell in real * space. */ virtual Tensor<1,dim> shape_grad (const unsigned int i, const Point &p) const; + /** + * Just like for @p{shape_grad}, + * but this function will be + * called when the shape function + * has more than one non-zero + * vector component. In that + * case, this function should + * return the gradient of the + * @p{component}-th vector + * component of the @p{i}th shape + * function at point @p{p}. + */ + virtual Tensor<1,dim> shape_grad_component (const unsigned int i, + const Point &p, + const unsigned int component) const; + /** * Return the tensor of second * derivatives of the @p{i}th @@ -455,18 +528,50 @@ class FiniteElementBase : public Subscriptor, * on the unit cell. The * derivatives are derivatives on * the unit cell with respect to - * unit cell coordinates. + * unit cell coordinates. If + * the finite element is + * vector-valued, then return the + * value of the only non-zero + * component of the vector value + * of this shape function. If the + * shape function has more than + * one non-zero component (which + * we refer to with the term + * non-primitive), then derived + * classes implementing this + * function should throw an + * exception of type + * @p{ExcShapeFunctionNotPrimitive}. In + * that case, use the + * @ref{shape_grad_grad_component} + * function. * * An * @p{ExcUnitShapeValuesDoNotExist} * is thrown if the shape values * of the @p{FiniteElement} under - * consideration depend on the + * consideration depends on the * shape of the cell in real * space. */ virtual Tensor<2,dim> shape_grad_grad (const unsigned int i, - const Point &p) const; + const Point &p) const; + + /** + * Just like for @p{shape_grad_grad}, + * but this function will be + * called when the shape function + * has more than one non-zero + * vector component. In that + * case, this function should + * return the gradient of the + * @p{component}-th vector + * component of the @p{i}th shape + * function at point @p{p}. + */ + virtual Tensor<2,dim> shape_grad_grad_component (const unsigned int i, + const Point &p, + const unsigned int component) const; /** * Projection from a fine grid @@ -553,81 +658,162 @@ class FiniteElementBase : public Subscriptor, const FullMatrix & prolongate (const unsigned int child) const; /** - * Return a readonly reference to the - * matrix which describes the constraints - * at the interface between a refined and - * an unrefined cell. + * Return a readonly reference to + * the matrix which describes the + * constraints at the interface + * between a refined and an + * unrefined cell. * - * The matrix is obviously empty in only - * one space dimension, since there are no - * constraints then. + * The matrix is obviously empty + * in only one space dimension, + * since there are no constraints + * then. */ const FullMatrix & constraints () const; /** - * Comparison operator. We also check for - * equality of the constraint matrix, - * which is quite an expensive operation. - * Do therefore - * use this function with care, if possible - * only for debugging purposes. + * Comparison operator. We also + * check for equality of the + * constraint matrix, which is + * quite an expensive operation. + * Do therefore use this function + * with care, if possible only + * for debugging purposes. * - * Since this function is not that important, - * we avoid an implementational question - * about comparing arrays and do not compare - * the matrix arrays @p{restriction} and + * Since this function is not + * that important, we avoid an + * implementational question + * about comparing arrays and do + * not compare the matrix arrays + * @p{restriction} and * @p{prolongation}. */ bool operator == (const FiniteElementBase &) const; /** - * Compute system index from components. + * Given a vector component and + * an index of a shape function + * within the shape functions + * corresponding to this vector + * component, return the index of + * this shape function within the + * shape functions of this + * element. If this is a scalar + * element, then the given + * component may only be zero, + * and the given component index + * is also the return value. + * + * If the finite element is + * vector-valued and has + * non-primitive shape functions, + * i.e. some of its shape + * functions are non-zero in more + * than just one vector + * component, then this function + * cannot be used since shape + * functions are no more + * associated with individual + * vector components, and an + * exception of type + * @p{ExcFENotPrimitive} is + * thrown. */ unsigned int component_to_system_index (const unsigned int component, const unsigned int component_index) const; /** - * Compute component and index from - * system index. - * - * Return value contains first - * component and second index in - * component. - */ - std::pair - system_to_component_index (const unsigned int index) const; - - /** - * Compute system index from components on a face. + * Same as above, but compute the + * data from the index of a shape + * function on a face. */ unsigned int face_component_to_system_index (const unsigned int component, const unsigned int component_index) const; - + /** - * Compute component and index from system - * index for a face. + * Compute vector component and + * index of this shape function + * within the shape functions + * corresponding to this + * component from the index of a + * shape function within this + * finite element. * - * Return value contains first - * component and second index in - * component. + * If the element is scalar, then + * the component is always zero, + * and the index within this + * component is equal to the + * overall index. + * + * If the shape function + * referenced has more than one + * non-zero component, then it + * cannot be associated with one + * vector component, and an + * exception of type + * @p{ExcShapeFunctionNotPrimitive} + * will be raised. + * + * Note that if the element is + * composed of other (base) + * elements, and a base element + * has more than one component + * but all its shape functions + * are primitive (i.e. are + * non-zero in only one + * component), then this mapping + * contains valid + * information. However, the + * index of a shape function of + * this element within one + * component (i.e. the second + * number of the respective entry + * of this array) does not + * indicate the index of the + * respective shape function + * within the base element (since + * that has more than one + * vector-component). For this + * information, refer to the + * @p{system_to_base_table} + * field. + */ + std::pair + system_to_component_index (const unsigned int index) const; + + /** + * Same as above, but do it for + * shape functions and their + * indices on a face. */ std::pair face_system_to_component_index (const unsigned int index) const; /** - * The base element establishing a - * component. + * Given a vector component, + * return an index which base + * element implements this + * component, and which vector + * component is this base element + * this is. This information is + * only of interest for + * vector-valued finite elements + * which are composed of several + * sub-elements. In that case, + * one may want to obtain + * information about the element + * implementing a certain vector + * component, which can be done + * using this function and the + * @ref{FESystem}::@p{base_element} + * function. * - * This table converts a - * component number to the - * @p{base_element} number. While - * component information contains - * multiplicity of base elements, - * the result allows access to - * shape functions of the base - * element. + * If this is a scalar finite + * element, then the return value + * is always equalt to zero. */ - unsigned int component_to_base (unsigned int index) const; + std::pair + component_to_base (unsigned int component) const; /** * Access the @p{restriction_is_additive_flag} @@ -746,6 +932,97 @@ class FiniteElementBase : public Subscriptor, * yields a non-empty array. */ bool has_face_support_points () const; + + /** + * Return in which of the vector + * components of this finite + * element the @p{i}the shape + * function is non-zero. The + * length of the returned array + * is equal to the number of + * vector components of this + * element. + * + * For most finite element + * spaces, the result of this + * function will be a vector with + * exactly one element being + * @p{true}, since for most + * spaces the individual vector + * components are + * independent. Only for those + * spaces that couple the + * components, for example to + * make a shape function + * divergence free, will there be + * more than one @p{true} entry. + */ + const std::vector & + get_nonzero_components (const unsigned int i) const; + + /** + * Return in how many vector + * components the @p{i}th shape + * function is non-zero. This + * value equals the number of + * entries equal to @p{true} in + * the result of the + * @p{get_nonzero_components} + * function. + * + * For most finite element + * spaces, the result will be + * equal to one. It is not equal + * to one only for those ansatz + * spaces for which vector-valued + * shape functions couple the + * individual components, for + * example in order to make them + * divergence-free. + */ + unsigned int + n_nonzero_components (const unsigned int i) const; + + /** + * Return whether the @p{i}th + * shape function is primitive in + * the sense that the shape + * function is non-zero in only + * one vector + * component. Non-primitive shape + * functions would then, for + * example, be those of + * divergence free ansatz spaces, + * in which the individual vector + * components are coupled. + * + * The result of the function is + * @p{true} if and only if the + * result of + * @p{n_nonzero_components(i)} is + * equal to one. + */ + bool + is_primitive (const unsigned int i) const; + + /** + * Return whether the entire + * finite element is primitive, + * in the sense that all its + * shape functions are + * primitive. If the finite + * element is scalar, then this + * is always the case. + * + * Since this is an extremely + * common operations, the result + * is cached in the + * @p{cached_primitivity} + * variable which is computed in + * the constructor. + */ + bool + is_primitive () const; /** * Determine an estimate for the @@ -760,6 +1037,21 @@ class FiniteElementBase : public Subscriptor, */ unsigned int memory_consumption () const; + /** + * Exception + */ + DeclException1 (ExcShapeFunctionNotPrimitive, + int, + << "The shape function with index " << arg1 + << " is not primitive, i.e. it is vector-valued and " + << "has more than one non-zero vector component. This " + << "function cannot be called for these shape functions. " + << "Maybe you want to use the same function with the " + << "_component suffix?"); + /** + * Exception + */ + DeclException0 (ExcFENotPrimitive); /** * Exception */ @@ -843,20 +1135,83 @@ class FiniteElementBase : public Subscriptor, FullMatrix interface_constraints; /** - * Map between linear dofs and - * component dofs. + * Store what + * @p{system_to_component_index} + * will return. */ std::vector< std::pair > system_to_component_table; /** - * Map between linear dofs and - * component dofs on face. + * Map between linear dofs and + * component dofs on face. This + * is filled with default values + * in the constructor, but + * derived classes will have to + * overwrite the information if + * necessary. + * + * By component, we mean the + * vector component, not the base + * element. The information thus + * makes only sense if a shape + * function is non-zero in only + * one component. */ std::vector< std::pair > face_system_to_component_table; + /** + * For each shape function, store + * to which base element and + * which instance of this base + * element (in case its + * multiplicity is greater than + * one) it belongs, and its index + * within this base element. If + * the element is not composed of + * others, then base and instance + * are always zero, and the index + * is equal to the number of the + * shape function. If the element + * is composed of single + * instances of other elements + * (i.e. all with multiplicity + * one) all of which are scalar, + * then base values and dof + * indices within this element + * are equal to the + * @p{system_to_component_table}. It + * differs only in case the + * element is composed of other + * elements and at least one of + * them is vector-valued itself. + * + * This array has valid values + * also in the case of + * vector-value + * (i.e. non-primitive) shape + * functions, in contrast to the + * @p{system_to_component_table}. + */ + std::vector,unsigned int> > system_to_base_table; + + /** + * Likewise for the indices on + * faces. + */ + std::vector,unsigned int> > face_system_to_base_table; + /** * Map between component and - * linear dofs. + * linear dofs: For each pair of + * vector component and index + * within this component, store + * the global dof number in the + * composed element. If the + * element is scalar, then the + * outer (component) index can + * only be zero, and the inner + * index is equal to the stored + * value. */ std::vector< std::vector > component_to_system_table; @@ -871,9 +1226,12 @@ class FiniteElementBase : public Subscriptor, * a component. * * This table converts a - * component number to the - * @p{base_element} number. While - * component information contains + * component number to a pair + * consisting of the + * @p{base_element} number, and + * the component within this base + * element. While component + * information contains * multiplicity of base elements, * the result allows access to * shape functions of the base @@ -890,8 +1248,8 @@ class FiniteElementBase : public Subscriptor, * case, the initialization by * the base class is sufficient. */ - std::vector component_to_base_table; - + std::vector > component_to_base_table; + /** * Projection matrices are * concatenated or summed up. @@ -957,6 +1315,48 @@ class FiniteElementBase : public Subscriptor, */ typename std::vector > unit_face_support_points; + /** + * For each shape function, give + * a vector of bools (with size + * equal to the number of vector + * components which this finite + * element has) indicating in + * which component each of these + * shape functions is non-zero. + * + * For primitive elements, there + * is only one non-zero + * component. + */ + const std::vector > nonzero_components; + + /** + * This array holds how many + * values in the respective entry + * of the @p{nonzero_components} + * element are non-zero. The + * array is thus a short-cut to + * allow faster access to this + * information than if we had to + * count the non-zero entries + * upon each request for this + * information. The field is + * initialized in the constructor + * of this class. + */ + const std::vector n_nonzero_components_table; + + /** + * Store whether all shape + * functions are primitive. Since + * finding this out is a very + * common operation, we cache the + * result, i.e. compute the value + * in the constructor for simpler + * access. + */ + const bool cached_primitivity; + /** * Compute second derivatives by * finite differences of @@ -982,7 +1382,20 @@ class FiniteElementBase : public Subscriptor, * 1e-6. */ static const double fd_step_length; - + + /** + * Given the pattern of nonzero + * components for each shape + * function, compute for each + * entry how many components are + * non-zero for each shape + * function. This function is + * used in the constructor of + * this class. + */ + static + std::vector + compute_n_nonzero_components (const std::vector > &nonzero_components); /** * Allow the FESystem class to @@ -1000,7 +1413,7 @@ class FiniteElementBase : public Subscriptor, * that throw if the matrices are * not already initialized. */ - friend class FESystem; + template friend class FESystem; }; @@ -1078,6 +1491,9 @@ FiniteElementBase::component_to_system_index (const unsigned int component, Assert(component_index::ExcFENotPrimitive()); + return component_to_system_table[component][component_index]; } @@ -1087,8 +1503,10 @@ inline std::pair FiniteElementBase::system_to_component_index (const unsigned int index) const { - Assert(index < system_to_component_table.size(), + Assert (index < system_to_component_table.size(), ExcIndexRange(index, 0, system_to_component_table.size())); + Assert (is_primitive (index), + typename FiniteElementBase::ExcShapeFunctionNotPrimitive(index)); return system_to_component_table[index]; } @@ -1104,6 +1522,9 @@ FiniteElementBase::face_component_to_system_index (const unsigned int compo Assert(component_index::ExcFENotPrimitive()); + return face_component_to_system_table[component][component_index]; } @@ -1115,19 +1536,21 @@ FiniteElementBase::face_system_to_component_index (const unsigned int index { Assert(index < face_system_to_component_table.size(), ExcIndexRange(index, 0, face_system_to_component_table.size())); +//TODO: check for primitivity of this shape function. this needs the global dof index +// Assert (is_primitive (face_to_cell_index(index)), +// typename FiniteElementBase::ExcShapeFunctionNotPrimitive(index)); return face_system_to_component_table[index]; } template inline -unsigned int +std::pair FiniteElementBase::component_to_base (unsigned int index) const { - if (n_components() == 1) - return 0; Assert(index < component_to_base_table.size(), ExcIndexRange(index, 0, component_to_base_table.size())); + return component_to_base_table[index]; } @@ -1143,4 +1566,57 @@ FiniteElementBase::restriction_is_additive (const unsigned int component) c } +template +inline +const std::vector & +FiniteElementBase::get_nonzero_components (const unsigned int i) const +{ + Assert (i < dofs_per_cell, ExcIndexRange (i, 0, dofs_per_cell)); + return nonzero_components[i]; +}; + + + +template +inline +unsigned int +FiniteElementBase::n_nonzero_components (const unsigned int i) const +{ + Assert (i < dofs_per_cell, ExcIndexRange (i, 0, dofs_per_cell)); + return n_nonzero_components_table[i]; +}; + + + +template +inline +bool +FiniteElementBase::is_primitive (const unsigned int i) const +{ + Assert (i < this->dofs_per_cell, ExcIndexRange (i, 0, dofs_per_cell)); + + // return primitivity of a shape + // function by checking whether it + // has more than one non-zero + // component or not. we could cache + // this value in an array of bools, + // but accessing a bit-vector (as + // std::vector is) is + // probably more expensive than + // just comparing against 1 + return (n_nonzero_components_table[i] == 1); +}; + + +template +inline +bool +FiniteElementBase::is_primitive () const +{ + return cached_primitivity; +}; + + + + #endif diff --git a/deal.II/deal.II/include/fe/fe_dgp.h b/deal.II/deal.II/include/fe/fe_dgp.h index 4dc9e2af9f..d6bcd24230 100644 --- a/deal.II/deal.II/include/fe/fe_dgp.h +++ b/deal.II/deal.II/include/fe/fe_dgp.h @@ -15,6 +15,7 @@ #include #include +#include #include template class PolynomialSpace; @@ -39,46 +40,111 @@ class FE_DGP : public FiniteElement * Constructor for tensor product * polynomials of degree @p{k}. */ - FE_DGP (unsigned int k); - /** - * Destructor. - */ - ~FE_DGP (); + FE_DGP (const unsigned int k); /** * Return the value of the * @p{i}th shape function at the - * point @p{p}. @p{p} is a point - * on the reference element. + * point @p{p}. See the + * @ref{FiniteElementBase} base + * class for more information + * about the semantics of this + * function. */ virtual double shape_value (const unsigned int i, const Point &p) const; + + /** + * Return the value of the + * @p{component}th vector + * component of the @p{i}th shape + * function at the point + * @p{p}. See the + * @ref{FiniteElementBase} base + * class for more information + * about the semantics of this + * function. + * + * Since this element is scalar, + * the returned value is the same + * as if the function without the + * @p{_component} suffix were + * called, provided that the + * specified component is zero. + */ + virtual double shape_value_component (const unsigned int i, + const Point &p, + const unsigned int component) const; /** * Return the gradient of the * @p{i}th shape function at the - * point @p{p}. @p{p} is a point - * on the reference element, and - * likewise the gradient is the - * gradient on the unit cell with - * respect to unit cell - * coordinates. + * point @p{p}. See the + * @ref{FiniteElementBase} base + * class for more information + * about the semantics of this + * function. */ virtual Tensor<1,dim> shape_grad (const unsigned int i, const Point &p) const; + /** + * Return the gradient of the + * @p{component}th vector + * component of the @p{i}th shape + * function at the point + * @p{p}. See the + * @ref{FiniteElementBase} base + * class for more information + * about the semantics of this + * function. + * + * Since this element is scalar, + * the returned value is the same + * as if the function without the + * @p{_component} suffix were + * called, provided that the + * specified component is zero. + */ + virtual Tensor<1,dim> shape_grad_component (const unsigned int i, + const Point &p, + const unsigned int component) const; + /** * Return the tensor of second * derivatives of the @p{i}th * shape function at point @p{p} - * on the unit cell. The - * derivatives are derivatives on - * the unit cell with respect to - * unit cell coordinates. + * on the unit cell. See the + * @ref{FiniteElementBase} base + * class for more information + * about the semantics of this + * function. */ virtual Tensor<2,dim> shape_grad_grad (const unsigned int i, const Point &p) const; + /** + * Return the second derivative + * of the @p{component}th vector + * component of the @p{i}th shape + * function at the point + * @p{p}. See the + * @ref{FiniteElementBase} base + * class for more information + * about the semantics of this + * function. + * + * Since this element is scalar, + * the returned value is the same + * as if the function without the + * @p{_component} suffix were + * called, provided that the + * specified component is zero. + */ + virtual Tensor<2,dim> shape_grad_grad_component (const unsigned int i, + const Point &p, + const unsigned int component) const; + /** * Return the polynomial degree * of this finite element, @@ -201,12 +267,12 @@ class FE_DGP : public FiniteElement /** * Declare a nested class which - * will has static definitions of + * will hold static definitions of * various matrices such as * constraint and embedding * matrices. The definition of * the various static fields are - * in the files @p{fe_q_[123]d.cc} + * in the files @p{fe_dgp_[123]d.cc} * in the source directory. */ struct Matrices @@ -259,12 +325,41 @@ class FE_DGP : public FiniteElement static std::vector get_dpo_vector(unsigned int degree); /** - * Compute flags for initial update only. + * Given a set of flags indicating + * what quantities are requested + * from a @p{FEValues} object, + * return which of these can be + * precomputed once and for + * all. Often, the values of + * shape function at quadrature + * points can be precomputed, for + * example, in which case the + * return value of this function + * would be the logical and of + * the input @p{flags} and + * @p{update_values}. + * + * For the present kind of finite + * element, this is exactly the + * case. */ virtual UpdateFlags update_once (const UpdateFlags flags) const; /** - * Compute flags for update on each cell. + * This is the opposite to the + * above function: given a set of + * flags indicating what we want + * to know, return which of these + * need to be computed each time + * we visit a new cell. + * + * If for the computation of one + * quantity something else is + * also required (for example, we + * often need the covariant + * transformation when gradients + * need to be computed), include + * this in the result as well. */ virtual UpdateFlags update_each (const UpdateFlags flags) const; @@ -274,13 +369,19 @@ class FE_DGP : public FiniteElement const unsigned int degree; /** - * Pointer to the tensor - * product polynomials. + * Pointer to an object + * representing the polynomial + * space used here. */ - PolynomialSpace* poly; + const PolynomialSpace polynomial_space; /** * Fields of cell-independent data. + * + * For information about the + * general purpose of this class, + * see the documentation of the + * base class. */ class InternalData : public FiniteElementBase::InternalDataBase { @@ -293,6 +394,17 @@ class FE_DGP : public FiniteElement * function, containing * values for each quadrature * point. + * + * In this array, we store + * the values of the shape + * function in the quadrature + * points on the unit + * cell. Since these values + * do not change under + * transformation to the real + * cell, we only need to copy + * them over when visiting a + * concrete cell. */ std::vector > shape_values; @@ -304,6 +416,15 @@ class FE_DGP : public FiniteElement * function, containing * values for each quadrature * point. + * + * We store the gradients in + * the quadrature points on + * the unit cell. We then + * only have to apply the + * transformation (which is a + * matrix-vector + * multiplication) when + * visiting an actual cell. */ typename std::vector > > shape_gradients; }; diff --git a/deal.II/deal.II/include/fe/fe_dgq.h b/deal.II/deal.II/include/fe/fe_dgq.h index d452fa0276..131b0a92dd 100644 --- a/deal.II/deal.II/include/fe/fe_dgq.h +++ b/deal.II/deal.II/include/fe/fe_dgq.h @@ -15,6 +15,7 @@ #include #include +#include #include template class TensorProductPolynomials; @@ -39,46 +40,111 @@ class FE_DGQ : public FiniteElement * Constructor for tensor product * polynomials of degree @p{k}. */ - FE_DGQ (unsigned int k); - /** - * Destructor. - */ - ~FE_DGQ (); + FE_DGQ (const unsigned int k); /** * Return the value of the * @p{i}th shape function at the - * point @p{p}. @p{p} is a point - * on the reference element. + * point @p{p}. See the + * @ref{FiniteElementBase} base + * class for more information + * about the semantics of this + * function. */ virtual double shape_value (const unsigned int i, const Point &p) const; + /** + * Return the value of the + * @p{component}th vector + * component of the @p{i}th shape + * function at the point + * @p{p}. See the + * @ref{FiniteElementBase} base + * class for more information + * about the semantics of this + * function. + * + * Since this element is scalar, + * the returned value is the same + * as if the function without the + * @p{_component} suffix were + * called, provided that the + * specified component is zero. + */ + virtual double shape_value_component (const unsigned int i, + const Point &p, + const unsigned int component) const; + /** * Return the gradient of the * @p{i}th shape function at the - * point @p{p}. @p{p} is a point - * on the reference element, and - * likewise the gradient is the - * gradient on the unit cell with - * respect to unit cell - * coordinates. + * point @p{p}. See the + * @ref{FiniteElementBase} base + * class for more information + * about the semantics of this + * function. */ virtual Tensor<1,dim> shape_grad (const unsigned int i, const Point &p) const; + /** + * Return the gradient of the + * @p{component}th vector + * component of the @p{i}th shape + * function at the point + * @p{p}. See the + * @ref{FiniteElementBase} base + * class for more information + * about the semantics of this + * function. + * + * Since this element is scalar, + * the returned value is the same + * as if the function without the + * @p{_component} suffix were + * called, provided that the + * specified component is zero. + */ + virtual Tensor<1,dim> shape_grad_component (const unsigned int i, + const Point &p, + const unsigned int component) const; + /** * Return the tensor of second * derivatives of the @p{i}th * shape function at point @p{p} - * on the unit cell. The - * derivatives are derivatives on - * the unit cell with respect to - * unit cell coordinates. + * on the unit cell. See the + * @ref{FiniteElementBase} base + * class for more information + * about the semantics of this + * function. */ virtual Tensor<2,dim> shape_grad_grad (const unsigned int i, const Point &p) const; + /** + * Return the second derivative + * of the @p{component}th vector + * component of the @p{i}th shape + * function at the point + * @p{p}. See the + * @ref{FiniteElementBase} base + * class for more information + * about the semantics of this + * function. + * + * Since this element is scalar, + * the returned value is the same + * as if the function without the + * @p{_component} suffix were + * called, provided that the + * specified component is zero. + */ + virtual Tensor<2,dim> shape_grad_grad_component (const unsigned int i, + const Point &p, + const unsigned int component) const; + /** * Return the polynomial degree * of this finite element, @@ -201,12 +267,12 @@ class FE_DGQ : public FiniteElement /** * Declare a nested class which - * will has static definitions of + * will hold static definitions of * various matrices such as * constraint and embedding * matrices. The definition of * the various static fields are - * in the files @p{fe_q_[123]d.cc} + * in the files @p{fe_dgq_[123]d.cc} * in the source directory. */ struct Matrices @@ -259,12 +325,41 @@ class FE_DGQ : public FiniteElement static std::vector get_dpo_vector(unsigned int degree); /** - * Compute flags for initial update only. + * Given a set of flags indicating + * what quantities are requested + * from a @p{FEValues} object, + * return which of these can be + * precomputed once and for + * all. Often, the values of + * shape function at quadrature + * points can be precomputed, for + * example, in which case the + * return value of this function + * would be the logical and of + * the input @p{flags} and + * @p{update_values}. + * + * For the present kind of finite + * element, this is exactly the + * case. */ virtual UpdateFlags update_once (const UpdateFlags flags) const; /** - * Compute flags for update on each cell. + * This is the opposite to the + * above function: given a set of + * flags indicating what we want + * to know, return which of these + * need to be computed each time + * we visit a new cell. + * + * If for the computation of one + * quantity something else is + * also required (for example, we + * often need the covariant + * transformation when gradients + * need to be computed), include + * this in the result as well. */ virtual UpdateFlags update_each (const UpdateFlags flags) const; @@ -306,10 +401,15 @@ class FE_DGQ : public FiniteElement * Pointer to the tensor * product polynomials. */ - TensorProductPolynomials* poly; + const TensorProductPolynomials polynomial_space; /** * Fields of cell-independent data. + * + * For information about the + * general purpose of this class, + * see the documentation of the + * base class. */ class InternalData : public FiniteElementBase::InternalDataBase { @@ -322,6 +422,17 @@ class FE_DGQ : public FiniteElement * function, containing * values for each quadrature * point. + * + * In this array, we store + * the values of the shape + * function in the quadrature + * points on the unit + * cell. Since these values + * do not change under + * transformation to the real + * cell, we only need to copy + * them over when visiting a + * concrete cell. */ std::vector > shape_values; @@ -333,6 +444,15 @@ class FE_DGQ : public FiniteElement * function, containing * values for each quadrature * point. + * + * We store the gradients in + * the quadrature points on + * the unit cell. We then + * only have to apply the + * transformation (which is a + * matrix-vector + * multiplication) when + * visiting an actual cell. */ typename std::vector > > shape_gradients; }; @@ -347,7 +467,7 @@ class FE_DGQ : public FiniteElement * access to build_renumbering * function. */ - friend class MappingQ; + template friend class MappingQ; }; diff --git a/deal.II/deal.II/include/fe/fe_nedelec.h b/deal.II/deal.II/include/fe/fe_nedelec.h new file mode 100644 index 0000000000..d5e8f390e4 --- /dev/null +++ b/deal.II/deal.II/include/fe/fe_nedelec.h @@ -0,0 +1,520 @@ +//--------------------------------------------------------------- +// $Id$ +// Version: $Name$ +// +// Copyright (C) 2002 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. +// +//--------------------------------------------------------------- +#ifndef __deal2__fe_nedelec_h +#define __deal2__fe_nedelec_h + +#include +#include +#include +#include + +template class TensorProductPolynomials; +template class MappingQ; + + + +/** + * Implementation of continuous Nedelec elements for the space H_curl. + * + * The constructor of this class takes the degree @p{p} of this finite + * element. + * + * + * @sect3{Numbering of the degrees of freedom (DoFs)} + * + * Nedelec elements have their degrees of freedom on edges, with shape + * functions being vector valued and pointing in tangential + * direction. We use the standard enumeration and direction of edges + * in deal.II, yielding the following shape functions in 2d: + * + * @begin{verbatim} + * 2 + * *--->---* + * | | + * 3^ ^1 + * | | + * *--->---* + * 0 + * @end{verbatim} + * + * For the 3d case, the ordering follows the same scheme: the lines + * are numbered as described in the documentation of the + * @ref{Triangulation} class, i.e. + * @begin{verbatim} + * *---6---* *---6---* + * /| | / /| + * 11 | 5 11 10 5 + * / 7 | / / | + * * | | *---2---* | + * | *---4---* | | * + * | / / | 1 / + * 3 8 9 3 | 9 + * |/ / | |/ + * *---0---* *---0---* + * @end{verbatim} + * and their directions are as follows: + * @begin{verbatim} + * *--->---* *--->---* + * /| | / /| + * ^ | ^ ^ ^ ^ + * / ^ | / / | + * * | | *--->---* | + * | *--->---* | | * + * | / / | ^ / + * ^ ^ ^ ^ | ^ + * |/ / | |/ + * *--->---* *--->---* + * @end{verbatim} + * + * The element does not make much sense in 1d, so it is not + * implemented there. + * + * + * @author Anna Schneebeli, Wolfgang Bangerth, 2002 + */ +template +class FE_Nedelec : public FiniteElement +{ + public: + /** + * Constructor for the Nedelec + * element of degree @p{p}. + */ + FE_Nedelec (const unsigned int p); + + /** + * Return the value of the + * @p{component}th vector + * component of the @p{i}th shape + * function at the point + * @p{p}. See the + * @ref{FiniteElementBase} base + * class for more information + * about the semantics of this + * function. + */ + virtual double shape_value_component (const unsigned int i, + const Point &p, + const unsigned int component) const; + + /** + * Return the gradient of the + * @p{component}th vector + * component of the @p{i}th shape + * function at the point + * @p{p}. See the + * @ref{FiniteElementBase} base + * class for more information + * about the semantics of this + * function. + */ + virtual Tensor<1,dim> shape_grad_component (const unsigned int i, + const Point &p, + const unsigned int component) const; + + /** + * Return the second derivative + * of the @p{component}th vector + * component of the @p{i}th shape + * function at the point + * @p{p}. See the + * @ref{FiniteElementBase} base + * class for more information + * about the semantics of this + * function. + */ + virtual Tensor<2,dim> shape_grad_grad_component (const unsigned int i, + const Point &p, + const unsigned int component) const; + + /** + * Return the polynomial degree + * of this finite element, + * i.e. the value passed to the + * constructor. + */ + unsigned int get_degree () const; + + /** + * Number of base elements in a + * mixed discretization. Here, + * this is of course equal to + * one. + */ + virtual unsigned int n_base_elements () const; + + /** + * Access to base element + * objects. Since this element is + * atomic, @p{base_element(0)} is + * @p{this}, and all other + * indices throw an error. + */ + virtual const FiniteElement & base_element (const unsigned int index) const; + + /** + * Check for non-zero values on a face. + * + * This function returns + * @p{true}, if the shape + * function @p{shape_index} has + * non-zero values on the face + * @p{face_index}. + * + * Implementation of the + * interface in + * @ref{FiniteElement} + */ + virtual bool has_support_on_face (const unsigned int shape_index, + const unsigned int face_index) const; + + /** + * Determine an estimate for the + * memory consumption (in bytes) + * of this object. + * + * This function is made virtual, + * since finite element objects + * are usually accessed through + * pointers to their base class, + * rather than the class itself. + */ + virtual unsigned int memory_consumption () const; + + /** + * Exception + */ + DeclException0 (ExcNotUsefulInThisDimension); + + protected: + /** + * @p{clone} function instead of + * a copy constructor. + * + * This function is needed by the + * constructors of @p{FESystem}. + */ + virtual FiniteElement * clone() const; + + /** + * Prepare internal data + * structures and fill in values + * independent of the cell. + */ + virtual + typename Mapping::InternalDataBase * + get_data (const UpdateFlags, + const Mapping& mapping, + const Quadrature& quadrature) const ; + + /** + * Implementation of the same + * function in + * @ref{FiniteElement}. + */ + virtual void + fill_fe_values (const Mapping &mapping, + const typename DoFHandler::cell_iterator &cell, + const Quadrature &quadrature, + typename Mapping::InternalDataBase &mapping_internal, + typename Mapping::InternalDataBase &fe_internal, + FEValuesData& data) const; + + /** + * Implementation of the same + * function in + * @ref{FiniteElement}. + */ + virtual void + fill_fe_face_values (const Mapping &mapping, + const typename DoFHandler::cell_iterator &cell, + const unsigned int face_no, + const Quadrature &quadrature, + typename Mapping::InternalDataBase &mapping_internal, + typename Mapping::InternalDataBase &fe_internal, + FEValuesData& data) const ; + + /** + * Implementation of the same + * function in + * @ref{FiniteElement}. + */ + virtual void + fill_fe_subface_values (const Mapping &mapping, + const typename DoFHandler::cell_iterator &cell, + const unsigned int face_no, + const unsigned int sub_no, + const Quadrature &quadrature, + typename Mapping::InternalDataBase &mapping_internal, + typename Mapping::InternalDataBase &fe_internal, + FEValuesData& data) const ; + + private: + + /** + * Declare a nested class which + * will hold static definitions of + * various matrices such as + * constraint and embedding + * matrices. The definition of + * the various static fields are + * in the files @p{fe_nedelec_[23]d.cc} + * in the source directory. + */ + struct Matrices + { + /** + * Embedding matrices. For + * each element type (the + * first index) there are as + * many embedding matrices as + * there are children per + * cell. The first index + * starts with linear + * elements and goes up in + * polynomial degree. The + * array may grow in the + * future with the number of + * elements for which these + * matrices have been + * computed. If for some + * element, the matrices have + * not been computed then you + * may use the element + * nevertheless but can not + * access the respective + * fields. + */ + static const double * const + embedding[][GeometryInfo::children_per_cell]; + + /** + * Number of elements (first + * index) the above field + * has. Equals the highest + * polynomial degree for + * which the embedding + * matrices have been + * computed. + */ + static const unsigned int n_embedding_matrices; + + /** + * As the + * @p{embedding_matrices} + * field, but for the + * interface constraints. One + * for each element for which + * it has been computed. + */ + static const double * const constraint_matrices[]; + + /** + * Like + * @p{n_embedding_matrices}, + * but for the number of + * interface constraint + * matrices. + */ + static const unsigned int n_constraint_matrices; + }; + + /** + * Only for internal use. Its + * full name is + * @p{get_dofs_per_object_vector} + * function and it creates the + * @p{dofs_per_object} vector that is + * needed within the constructor to + * be passed to the constructor of + * @p{FiniteElementData}. + */ + static std::vector get_dpo_vector(const unsigned int degree); + + /** + * Initialize the + * @p{unit_support_points} field + * of the @ref{FiniteElementBase} + * class. Called from the + * constructor. + */ + void initialize_unit_support_points (); + + /** + * Initialize the + * @p{unit_face_support_points} field + * of the @ref{FiniteElementBase} + * class. Called from the + * constructor. + */ + void initialize_unit_face_support_points (); + + /** + * Given a set of flags indicating + * what quantities are requested + * from a @p{FEValues} object, + * return which of these can be + * precomputed once and for + * all. Often, the values of + * shape function at quadrature + * points can be precomputed, for + * example, in which case the + * return value of this function + * would be the logical and of + * the input @p{flags} and + * @p{update_values}. + * + * For the present kind of finite + * element, this is exactly the + * case. + */ + virtual UpdateFlags update_once (const UpdateFlags flags) const; + + /** + * This is the opposite to the + * above function: given a set of + * flags indicating what we want + * to know, return which of these + * need to be computed each time + * we visit a new cell. + * + * If for the computation of one + * quantity something else is + * also required (for example, we + * often need the covariant + * transformation when gradients + * need to be computed), include + * this in the result as well. + */ + virtual UpdateFlags update_each (const UpdateFlags flags) const; + + /** + * Degree of the polynomials. + */ + const unsigned int degree; + + /** + * Fields of cell-independent data. + * + * For information about the + * general purpose of this class, + * see the documentation of the + * base class. + */ + class InternalData : public FiniteElementBase::InternalDataBase + { + public: + /** + * Array with shape function + * values in quadrature + * points. There is one + * vector for each shape + * function, containing + * values for each quadrature + * point. Since the shape + * functions are + * vector-valued (with as + * many components as there + * are space dimensions), the + * value is a tensor. + * + * In this array, we store + * the values of the shape + * function in the quadrature + * points on the unit + * cell. The transformation + * to the real space cell is + * then simply done by + * multiplication with the + * Jacobian of the mapping. + */ + std::vector > > shape_values; + + /** + * Array with shape function + * gradients in quadrature + * points. There is one + * vector for each shape + * function, containing + * values for each quadrature + * point. + * + * We store the gradients in + * the quadrature points on + * the unit cell. We then + * only have to apply the + * transformation (which is a + * matrix-vector + * multiplication) when + * visiting an actual cell. + */ + typename std::vector > > shape_gradients; + }; + + /** + * Allow access from other + * dimensions. + */ + template friend class FE_Nedelec; +}; + + +/* -------------- declaration of explicit specializations ------------- */ + +template <> void FE_Nedelec<1>::initialize_unit_face_support_points (); + +// declaration of explicit specializations of member variables, if the +// compiler allows us to do that (the standard says we must) +#ifndef DEAL_II_MEMBER_VAR_SPECIALIZATION_BUG +template <> +const double * const +FE_Nedelec<1>::Matrices::embedding[][GeometryInfo<1>::children_per_cell]; + +template <> +const unsigned int FE_Nedelec<1>::Matrices::n_embedding_matrices; + +template <> +const double * const FE_Nedelec<1>::Matrices::constraint_matrices[]; + +template <> +const unsigned int FE_Nedelec<1>::Matrices::n_constraint_matrices; + +template <> +const double * const +FE_Nedelec<2>::Matrices::embedding[][GeometryInfo<2>::children_per_cell]; + +template <> +const unsigned int FE_Nedelec<2>::Matrices::n_embedding_matrices; + +template <> +const double * const FE_Nedelec<2>::Matrices::constraint_matrices[]; + +template <> +const unsigned int FE_Nedelec<2>::Matrices::n_constraint_matrices; + +template <> +const double * const +FE_Nedelec<3>::Matrices::embedding[][GeometryInfo<3>::children_per_cell]; + +template <> +const unsigned int FE_Nedelec<3>::Matrices::n_embedding_matrices; + +template <> +const double * const FE_Nedelec<3>::Matrices::constraint_matrices[]; + +template <> +const unsigned int FE_Nedelec<3>::Matrices::n_constraint_matrices; + +#endif + +#endif diff --git a/deal.II/deal.II/include/fe/fe_q.h b/deal.II/deal.II/include/fe/fe_q.h index d8b9a85d25..61abd087d4 100644 --- a/deal.II/deal.II/include/fe/fe_q.h +++ b/deal.II/deal.II/include/fe/fe_q.h @@ -15,6 +15,7 @@ #include #include +#include #include template class TensorProductPolynomials; @@ -245,45 +246,109 @@ class FE_Q : public FiniteElement */ FE_Q (const unsigned int p); - /** - * Destructor. - */ - ~FE_Q (); - /** * Return the value of the * @p{i}th shape function at the - * point @p{p}. @p{p} is a point - * on the reference element. + * point @p{p}. See the + * @ref{FiniteElementBase} base + * class for more information + * about the semantics of this + * function. */ virtual double shape_value (const unsigned int i, const Point &p) const; + /** + * Return the value of the + * @p{component}th vector + * component of the @p{i}th shape + * function at the point + * @p{p}. See the + * @ref{FiniteElementBase} base + * class for more information + * about the semantics of this + * function. + * + * Since this element is scalar, + * the returned value is the same + * as if the function without the + * @p{_component} suffix were + * called, provided that the + * specified component is zero. + */ + virtual double shape_value_component (const unsigned int i, + const Point &p, + const unsigned int component) const; + /** * Return the gradient of the * @p{i}th shape function at the - * point @p{p}. @p{p} is a point - * on the reference element, and - * likewise the gradient is the - * gradient on the unit cell with - * respect to unit cell - * coordinates. + * point @p{p}. See the + * @ref{FiniteElementBase} base + * class for more information + * about the semantics of this + * function. */ virtual Tensor<1,dim> shape_grad (const unsigned int i, const Point &p) const; + /** + * Return the gradient of the + * @p{component}th vector + * component of the @p{i}th shape + * function at the point + * @p{p}. See the + * @ref{FiniteElementBase} base + * class for more information + * about the semantics of this + * function. + * + * Since this element is scalar, + * the returned value is the same + * as if the function without the + * @p{_component} suffix were + * called, provided that the + * specified component is zero. + */ + virtual Tensor<1,dim> shape_grad_component (const unsigned int i, + const Point &p, + const unsigned int component) const; + /** * Return the tensor of second * derivatives of the @p{i}th * shape function at point @p{p} - * on the unit cell. The - * derivatives are derivatives on - * the unit cell with respect to - * unit cell coordinates. + * on the unit cell. See the + * @ref{FiniteElementBase} base + * class for more information + * about the semantics of this + * function. */ virtual Tensor<2,dim> shape_grad_grad (const unsigned int i, const Point &p) const; + /** + * Return the second derivative + * of the @p{component}th vector + * component of the @p{i}th shape + * function at the point + * @p{p}. See the + * @ref{FiniteElementBase} base + * class for more information + * about the semantics of this + * function. + * + * Since this element is scalar, + * the returned value is the same + * as if the function without the + * @p{_component} suffix were + * called, provided that the + * specified component is zero. + */ + virtual Tensor<2,dim> shape_grad_grad_component (const unsigned int i, + const Point &p, + const unsigned int component) const; + /** * Return the polynomial degree * of this finite element, @@ -405,7 +470,7 @@ class FE_Q : public FiniteElement /** * Declare a nested class which - * will has static definitions of + * will hold static definitions of * various matrices such as * constraint and embedding * matrices. The definition of @@ -584,14 +649,41 @@ class FE_Q : public FiniteElement void initialize_unit_face_support_points (); /** - * Compute flags for initial - * update only. + * Given a set of flags indicating + * what quantities are requested + * from a @p{FEValues} object, + * return which of these can be + * precomputed once and for + * all. Often, the values of + * shape function at quadrature + * points can be precomputed, for + * example, in which case the + * return value of this function + * would be the logical and of + * the input @p{flags} and + * @p{update_values}. + * + * For the present kind of finite + * element, this is exactly the + * case. */ virtual UpdateFlags update_once (const UpdateFlags flags) const; /** - * Compute flags for update on - * each cell. + * This is the opposite to the + * above function: given a set of + * flags indicating what we want + * to know, return which of these + * need to be computed each time + * we visit a new cell. + * + * If for the computation of one + * quantity something else is + * also required (for example, we + * often need the covariant + * transformation when gradients + * need to be computed), include + * this in the result as well. */ virtual UpdateFlags update_each (const UpdateFlags flags) const; @@ -624,33 +716,70 @@ class FE_Q : public FiniteElement * Pointer to the tensor * product polynomials. */ - TensorProductPolynomials* poly; + const TensorProductPolynomials polynomial_space; /** * Fields of cell-independent data. + * + * For information about the + * general purpose of this class, + * see the documentation of the + * base class. */ class InternalData : public FiniteElementBase::InternalDataBase { public: /** - * Array with shape function values - * in quadrature points. There is one - * vector for each shape function, containing - * values for each quadrature point. + * Array with shape function + * values in quadrature + * points. There is one + * vector for each shape + * function, containing + * values for each quadrature + * point. + * + * In this array, we store + * the values of the shape + * function in the quadrature + * points on the unit + * cell. Since these values + * do not change under + * transformation to the real + * cell, we only need to copy + * them over when visiting a + * concrete cell. */ std::vector > shape_values; /** - * Array with shape function gradients - * in quadrature points. There is one - * vector for each shape function, containing - * values for each quadrature point. - */ + * Array with shape function + * gradients in quadrature + * points. There is one + * vector for each shape + * function, containing + * values for each quadrature + * point. + * + * We store the gradients in + * the quadrature points on + * the unit cell. We then + * only have to apply the + * transformation (which is a + * matrix-vector + * multiplication) when + * visiting an actual cell. + */ typename std::vector > > shape_gradients; }; /** - * Allow access from other dimensions. + * Allow access from other + * dimensions. We need this since + * we want to call the functions + * @p{get_dpo_vector} and + * @p{lexicographic_to_hierarchic_numbering} + * for the faces of the finite + * element of dimension dim+1. */ template friend class FE_Q; }; diff --git a/deal.II/deal.II/include/fe/fe_system.h b/deal.II/deal.II/include/fe/fe_system.h index ef94f30669..663c8a925e 100644 --- a/deal.II/deal.II/include/fe/fe_system.h +++ b/deal.II/deal.II/include/fe/fe_system.h @@ -28,7 +28,7 @@ * into one. To the outside world, the resulting object looks just * like a usual finite element object, which is composed of several * other finite elements that are possibly of different type. - + * * The overall numbering of degrees of freedom is as follows: for each * subobject (vertex, line, quad, or hex), the degrees of freedom are * numbered such that we run over all subelements first, before @@ -127,7 +127,18 @@ class FESystem : public FiniteElement * Return the value of the * @p{i}th shape function at the * point @p{p}. @p{p} is a point - * on the reference element. + * on the reference element. Since + * this finite element is always + * vector-valued, we return the + * value of the only non-zero + * component of the vector value + * of this shape function. If the + * shape function has more than + * one non-zero component (which + * we refer to with the term + * non-primitive), then throw an + * exception of type + * @p{ExcShapeFunctionNotPrimitive}. * * An * @p{ExcUnitShapeValuesDoNotExist} @@ -140,6 +151,26 @@ class FESystem : public FiniteElement */ virtual double shape_value (const unsigned int i, const Point &p) const; + + /** + * Return the value of the + * @p{component}th vector + * component of the @p{i}th shape + * function at the point + * @p{p}. See the + * @ref{FiniteElementBase} base + * class for more information + * about the semantics of this + * function. + * + * Since this element is vector + * valued in general, it relays + * the computation of these + * values to the base elements. + */ + virtual double shape_value_component (const unsigned int i, + const Point &p, + const unsigned int component) const; /** * Return the gradient of the @@ -149,7 +180,18 @@ class FESystem : public FiniteElement * likewise the gradient is the * gradient on the unit cell with * respect to unit cell - * coordinates. + * coordinates. Since + * this finite element is always + * vector-valued, we return the + * value of the only non-zero + * component of the vector value + * of this shape function. If the + * shape function has more than + * one non-zero component (which + * we refer to with the term + * non-primitive), then throw an + * exception of type + * @p{ExcShapeFunctionNotPrimitive}. * * An * @p{ExcUnitShapeValuesDoNotExist} @@ -163,6 +205,26 @@ class FESystem : public FiniteElement virtual Tensor<1,dim> shape_grad (const unsigned int i, const Point &p) const; + /** + * Return the gradient of the + * @p{component}th vector + * component of the @p{i}th shape + * function at the point + * @p{p}. See the + * @ref{FiniteElementBase} base + * class for more information + * about the semantics of this + * function. + * + * Since this element is vector + * valued in general, it relays + * the computation of these + * values to the base elements. + */ + virtual Tensor<1,dim> shape_grad_component (const unsigned int i, + const Point &p, + const unsigned int component) const; + /** * Return the tensor of second * derivatives of the @p{i}th @@ -170,7 +232,18 @@ class FESystem : public FiniteElement * on the unit cell. The * derivatives are derivatives on * the unit cell with respect to - * unit cell coordinates. + * unit cell coordinates. Since + * this finite element is always + * vector-valued, we return the + * value of the only non-zero + * component of the vector value + * of this shape function. If the + * shape function has more than + * one non-zero component (which + * we refer to with the term + * non-primitive), then throw an + * exception of type + * @p{ExcShapeFunctionNotPrimitive}. * * An * @p{ExcUnitShapeValuesDoNotExist} @@ -184,6 +257,26 @@ class FESystem : public FiniteElement virtual Tensor<2,dim> shape_grad_grad (const unsigned int i, const Point &p) const; + /** + * Return the second derivatives + * of the @p{component}th vector + * component of the @p{i}th shape + * function at the point + * @p{p}. See the + * @ref{FiniteElementBase} base + * class for more information + * about the semantics of this + * function. + * + * Since this element is vector + * valued in general, it relays + * the computation of these + * values to the base elements. + */ + virtual Tensor<2,dim> shape_grad_grad_component (const unsigned int i, + const Point &p, + const unsigned int component) const; + /** * Number of different base * elements of this object. @@ -215,7 +308,8 @@ class FESystem : public FiniteElement * multiplicities are greater * than one. */ - virtual const FiniteElement & base_element (const unsigned int index) const; + virtual const FiniteElement & + base_element (const unsigned int index) const; /** * Check for non-zero values on a face. @@ -463,6 +557,51 @@ class FESystem : public FiniteElement const unsigned int N2, const FiniteElement &fe3, const unsigned int N3); + + /** + * Compute the non-zero vector + * components of a composed + * finite element. + */ + static std::vector > + compute_nonzero_components (const FiniteElement &fe1, + const unsigned int N1); + + /** + * Compute the non-zero vector + * components of a composed + * finite element. + */ + static std::vector > + compute_nonzero_components (const FiniteElement &fe1, + const unsigned int N1, + const FiniteElement &fe2, + const unsigned int N2); + + /** + * Compute the non-zero vector + * components of a composed + * finite element. + */ + static std::vector > + compute_nonzero_components (const FiniteElement &fe1, + const unsigned int N1, + const FiniteElement &fe2, + const unsigned int N2, + const FiniteElement &fe3, + const unsigned int N3); + + /** + * Compute the nonzero components + * of a list of finite elements + * with multiplicities given in + * the second argument. This + * function is called from all + * the above functions. + */ + static std::vector > + compute_nonzero_components (const std::vector*> &fes, + const std::vector &multiplicities); /** * This function is simply @@ -479,12 +618,12 @@ class FESystem : public FiniteElement /** * Used by @p{initialize}. */ - void build_cell_table(); + void build_cell_tables(); /** * Used by @p{initialize}. */ - void build_face_table(); + void build_face_tables(); /** * Used by @p{initialize}. @@ -495,12 +634,11 @@ class FESystem : public FiniteElement * Usually: Fields of * cell-independent data. * - * But for @p{FESystem} this - * @p{InternalData} class does + * However, here, this class does * not itself store the data but * only pointers to - * @p{InternalDatas} of the base - * elements. + * @p{InternalData} objects for + * each of the base elements. */ class InternalData : public FiniteElementBase::InternalDataBase { @@ -513,7 +651,7 @@ class FESystem : public FiniteElement * vector to * @p{n_base_elements}. */ - InternalData(const unsigned int n_base_elements); + InternalData (const unsigned int n_base_elements); /** * Destructor. Deletes all @@ -593,26 +731,33 @@ class FESystem : public FiniteElement private: /** - * Pointers to the - * @p{InternalDatas} of the - * base elements. They are + * Pointers to + * @p{InternalData} objects + * for each of the base + * elements. They are * accessed to by the * @p{set_} and - * @p{get_fe_data} - * functions. + * @p{get_fe_data} functions. * * The size of this vector is * set to @p{n_base_elements} * by the InternalData - * constructor. It is - * filled by the @p{get_data} - * function. + * constructor. It is filled + * by the @p{get_data} + * function. Note that since + * the data for each instance + * of a base class is + * necessarily the same, we + * only need as many of these + * objects as there are base + * elements, irrespective of + * their multiplicity. */ typename std::vector::InternalDataBase *> base_fe_datas; /** * Pointers to the - * @p{FEValuesDatas} + * @p{FEValuesData} objects * that are given to the * @p{fill_fe_values} * function of the base diff --git a/deal.II/deal.II/include/fe/fe_values.h b/deal.II/deal.II/include/fe/fe_values.h index 1f29f44e28..3da628e874 100644 --- a/deal.II/deal.II/include/fe/fe_values.h +++ b/deal.II/deal.II/include/fe/fe_values.h @@ -28,6 +28,8 @@ #include #include +#include + template class Quadrature; @@ -46,57 +48,112 @@ class FEValuesData { public: /** - * Initialize all vectors to correct size. + * Initialize all vectors to + * correct size. */ - void initialize (const unsigned int n_quadrature_points, - const unsigned int n_shapes, - const UpdateFlags flags); + void initialize (const unsigned int n_quadrature_points, + const FiniteElement &fe, + const UpdateFlags flags); /** - * Storage type for shape values. + * Storage type for shape + * values. Each row in the matrix + * denotes the values of a single + * shape function at the + * different points, columns are + * for a single point with the + * different shape functions. + * + * If a shape function has more + * than one non-zero component + * (in deal.II diction: it is + * non-primitive), then we + * allocate one row per non-zero + * component, and shift + * subsequent rows backward. + * Lookup of the correct row for + * a shape function is thus + * simple in case the entire + * finite element is primitive + * (i.e. all shape functions are + * primitive), since then the + * shape function number equals + * the row number. Otherwise, use + * the + * @p{shape_function_to_row_table} + * array to get at the first row + * that belongs to this + * particular shape function, and + * navigate among all the rows + * for this shape function using + * the + * @p{FiniteElement::get_nonzero_components} + * function which tells us which + * components are non-zero and + * thus have a row in the array + * presently under discussion. */ typedef FullMatrix ShapeVector; /** - * Storage type for gradients. + * Storage type for + * gradients. The layout of data + * is the same as for the + * @ref{ShapeVector} data type. */ - typedef typename std::vector > > GradientVector; + typedef + typename std::vector > > + GradientVector; + + /** + * Likewise for second order + * derivatives. + */ + typedef + typename std::vector > > + GradGradVector; /** - * Store the values of the shape functions - * at the quadrature points. Rows in the - * matrices denote the values of a single - * shape function at the different points, - * columns are for a single point with the - * different shape functions. + * Store the values of the shape + * functions at the quadrature + * points. See the description of + * the data type for the layout + * of the data in this field. */ ShapeVector shape_values; /** - * Store the gradients of the shape - * functions at the quadrature points. - * For each shape function, there is a vector - * of the gradients in each quadrature point. + * Store the gradients of the + * shape functions at the + * quadrature points. See the + * description of the data type + * for the layout of the data in + * this field. */ GradientVector shape_gradients; /** - * Store the 2nd derivatives of the shape - * functions at the quadrature points. - * For each shape function, there is a vector - * of the 2nd gradients in each quadrature point. + * Store the 2nd derivatives of + * the shape functions at the + * quadrature points. See the + * description of the data type + * for the layout of the data in + * this field. */ - typename std::vector > > shape_2nd_derivatives; + GradGradVector shape_2nd_derivatives; /** - * Store an array of weights times the - * Jacobi determinant at the quadrature - * points. This function is reset each time - * @p{reinit} is called. The Jacobi determinant - * is actually the reciprocal value of the - * Jacobi matrices stored in this class, - * see the general documentation of this - * class for more information. + * Store an array of weights + * times the Jacobi determinant + * at the quadrature points. This + * function is reset each time + * @p{reinit} is called. The + * Jacobi determinant is actually + * the reciprocal value of the + * Jacobi matrices stored in this + * class, see the general + * documentation of this class + * for more information. */ std::vector JxW_values; @@ -123,6 +180,36 @@ class FEValuesData */ typename std::vector > boundary_forms; + /** + * Indicate the first row which a + * given shape function occupies + * in the @p{ShapeVector}, + * @p{ShapeGradient}, etc + * arrays. If all shape functions + * are primitive, then this is + * the identity mapping. If, on + * the other hand some shape + * functions have more than one + * non-zero vector components, + * then they may occupy more than + * one row, and this array + * indicates which is the first + * one. + * + * The questions which particular + * vector component occupies + * which row for a given shape + * function is answered as + * follows: we allocate one row + * for each non-zero component as + * indicated by the + * @p{FiniteElement::get_nonzero_components()} + * function, and the rows are in + * ascending order exactly those + * non-zero components. + */ + std::vector shape_function_to_row_table; + /** * Original update flags handed * to the constructor of @@ -151,7 +238,7 @@ class FEValuesData * called by the constructor and @p{reinit} functions of * @p{FEValues*}, respectively. * - * \subsection{General usage} + * @sect3{General usage} * * Usually, an object of @p{FEValues*} is used in integration loops * over all cells of a triangulation. To take full advantage of the @@ -160,7 +247,7 @@ class FEValuesData * like a magnifying glass being used to look at one item after the * other. A typical piece of code looks like this: * - * \begin{verbatim} + * @begin{verbatim} * FEValues values (mapping, finite_element, quadrature, flags); * for (cell = dof_handler.begin_active(); * cell != dof_handler.end(); @@ -169,30 +256,45 @@ class FEValuesData * values.reinit(cell); * ... * } - * \end{verbatim} + * @end{verbatim} * * - * \subsection{Member functions} + * @sect3{Member functions} * * The functions of this class fall into different cathegories: - * \begin{itemize} - * \item @p{shape_value}, @p{shape_grad}, etc: return one of the values + * @begin{itemize} + * @item @p{shape_value}, @p{shape_grad}, etc: return one of the values * of this object at a time. These functions are inlined, so this - * is the suggested access to all finite element values. There should be - * no loss in performance with an optimizing compiler. + * is the suggested access to all finite element values. There + * should be no loss in performance with an optimizing compiler. If + * the finite element is vector valued, then these functions return + * the only non-zero component of the requested shape + * function. However, some finite elements have shape functions + * that have more than one non-zero component (we call them + * non-"primitive"), and in this case this set of functions will + * throw an exception since they cannot generate a useful + * result. Rather, use the next set of functions. + * + * @item @p{shape_value_component}, @p{shape_grad_component}, etc: + * This is the same set of functions as above, except that for vector + * valued finite elements they return only one vector component. This + * is useful for elements of which shape functions have more than one + * non-zero component, since then the above functions cannot be used, + * and you have to walk over all (or only the non-zero) components of + * the shape function using this set of functions. * - * \item @p{get_function_values}, @p{get_function_grads}, @p{...}: + * @item @p{get_function_values}, @p{get_function_grads}, @p{...}: * Compute a finite element function or its derivative * in quadrature points. * - * \item @p{reinit}: initialize the @p{FEValues} object for a certain cell. + * @item @p{reinit}: initialize the @p{FEValues} object for a certain cell. * This function is not in the present class but only in the derived * classes and has a variable call syntax. * See the docs for the derived classes for more information. - * \end{itemize} + * @end{itemize} * * - * \subsection{UpdateFlags} + * @sect3{UpdateFlags} * * The @ref{UpdateFlags} object handed to the constructor is used to * determine, which of the data fields to compute. This way, it is @@ -217,12 +319,13 @@ class FEValuesBase : protected FEValuesData const unsigned int n_quadrature_points; /** - * Number of shape functions - * per cell. If we use this base class - * to evaluate a finite element on - * faces of cells, this is still the - * number of degrees of freedom per - * cell, not per face. + * Number of shape functions per + * cell. If we use this base + * class to evaluate a finite + * element on faces of cells, + * this is still the number of + * degrees of freedom per cell, + * not per face. */ const unsigned int dofs_per_cell; @@ -254,34 +357,57 @@ class FEValuesBase : protected FEValuesData ~FEValuesBase (); /** - * Value of the @p{function_no}th shape - * function at the @p{point_no}th quadrature + * Value of the @p{function_no}th + * shape function at the + * @p{point_no}th quadrature * point on the cell, face or * subface selected the last time * the @p{reinit} function of the * derived class was called. + * + * If the shape function is + * vector-valued, then this + * returns the only non-zero + * component. If the shape + * function has more than one + * non-zero component (i.e. it is + * not primitive), then throw an + * exception of type + * @p{ExcShapeFunctionNotPrimitive}. In + * that case, use the + * @ref{shape_value_component} + * function. */ double shape_value (const unsigned int function_no, const unsigned int point_no) const; /** - * This function is now - * deprecated and will be removed - * in the next release of the - * deal.II library. - * - * Pointer to the matrix holding - * all values of shape functions - * at all integration points, on - * the present cell, face or - * subface selected the last time - * the @p{reinit} function of the - * derived class was called. For - * the format of this matrix, see - * the documentation for the - * matrix itself. - */ - const typename FEValuesData::ShapeVector & get_shape_values () const; + * Return one vector component of + * the value of a shape function + * at a quadrature point. If the + * finite element is scalar, then + * only component zero is allowed + * and the return value equals + * that of the @p{shape_value} + * function. If the finite + * element is vector valued but + * all shape functions are + * primitive (i.e. they are + * non-zero in only one + * component), then the value + * returned by @p{shape_value} + * equals that of this function + * for exactly one + * component. This function is + * therefore only of greater + * interest if the shape function + * is not primitive, but then it + * is necessary since the other + * function cannot be used. + */ + double shape_value_component (const unsigned int function_no, + const unsigned int point_no, + const unsigned int component) const; /** * Values of the finite @@ -351,25 +477,52 @@ class FEValuesBase : protected FEValuesData * gradient's value is returned, * there should be no major * performance drawback. - */ - const Tensor<1,dim> & shape_grad (const unsigned int function, - const unsigned int quadrature_point) const; - - /** - * This function is now - * deprecated and will be removed - * in the next release of the - * deal.II library. * - * Pointer to the matrix holding - * all gradients of shape - * functions at all integration - * points, on the present cell. - * For the format of this matrix, - * see the documentation for the - * matrix itself. - */ - const typename FEValuesData::GradientVector & get_shape_grads () const; + * If the shape function is + * vector-valued, then this + * returns the only non-zero + * component. If the shape + * function has more than one + * non-zero component (i.e. it is + * not primitive), then throw an + * exception of type + * @p{ExcShapeFunctionNotPrimitive}. In + * that case, use the + * @ref{shape_grad_component} + * function. + */ + const Tensor<1,dim> & + shape_grad (const unsigned int function, + const unsigned int quadrature_point) const; + + /** + * Return one vector component of + * the gradient of a shape function + * at a quadrature point. If the + * finite element is scalar, then + * only component zero is allowed + * and the return value equals + * that of the @p{shape_grad} + * function. If the finite + * element is vector valued but + * all shape functions are + * primitive (i.e. they are + * non-zero in only one + * component), then the value + * returned by @p{shape_grad} + * equals that of this function + * for exactly one + * component. This function is + * therefore only of greater + * interest if the shape function + * is not primitive, but then it + * is necessary since the other + * function cannot be used. + */ + Tensor<1,dim> + shape_grad_component (const unsigned int function_no, + const unsigned int point_no, + const unsigned int component) const; /** * Gradients of the finite @@ -438,26 +591,56 @@ class FEValuesBase : protected FEValuesData * only a reference to the * derivative values is returned, * there should be no major - * performance drawback. */ - const Tensor<2,dim> & shape_2nd_derivative (const unsigned int function_no, - const unsigned int point_no) const; - - /** - * This function is now - * deprecated and will be removed - * in the next release of the - * deal.II library. + * performance drawback. * - * Pointer to the - * matrix holding all 2nd - * derivatives of shape functions - * at all integration points, on - * the present cell. For the - * format of this matrix, see the - * documentation for the matrix - * itself. - */ - const typename std::vector > > & get_shape_2nd_derivatives () const; + * If the shape function is + * vector-valued, then this + * returns the only non-zero + * component. If the shape + * function has more than one + * non-zero component (i.e. it is + * not primitive), then throw an + * exception of type + * @p{ExcShapeFunctionNotPrimitive}. In + * that case, use the + * @ref{shape_grad_grad_component} + * function. + */ + const Tensor<2,dim> & + shape_2nd_derivative (const unsigned int function_no, + const unsigned int point_no) const; + + + /** + * Return one vector component of + * the gradient of a shape + * function at a quadrature + * point. If the finite element + * is scalar, then only component + * zero is allowed and the return + * value equals that of the + * @p{shape_2nd_derivative} + * function. If the finite + * element is vector valued but + * all shape functions are + * primitive (i.e. they are + * non-zero in only one + * component), then the value + * returned by + * @p{shape_2nd_derivative} + * equals that of this function + * for exactly one + * component. This function is + * therefore only of greater + * interest if the shape function + * is not primitive, but then it + * is necessary since the other + * function cannot be used. + */ + Tensor<2,dim> + shape_2nd_derivative_component (const unsigned int function_no, + const unsigned int point_no, + const unsigned int component) const; /** * Tensor of second derivatives @@ -591,6 +774,21 @@ class FEValuesBase : protected FEValuesData * Exception */ DeclException0 (ExcFEDontMatch); + /** + * Exception + */ + DeclException1 (ExcShapeFunctionNotPrimitive, + int, + << "The shape function with index " << arg1 + << " is not primitive, i.e. it is vector-valued and " + << "has more than one non-zero vector component. This " + << "function cannot be called for these shape functions. " + << "Maybe you want to use the same function with the " + << "_component suffix?"); + /** + * Exception + */ + DeclException0 (ExcFENotPrimitive); protected: /** @@ -980,41 +1178,283 @@ double FEValuesBase::shape_value (const unsigned int i, const unsigned int j) const { - Assert (update_flags & update_values, ExcAccessToUninitializedField()); - return shape_values(i,j); + Assert (this->update_flags & update_values, + ExcAccessToUninitializedField()); + Assert (fe->is_primitive (i), + ExcShapeFunctionNotPrimitive(i)); + + // if the entire FE is primitive, + // then we can take a short-cut: + if (fe->is_primitive()) + return this->shape_values(i,j); + else + // otherwise, use the mapping + // between shape function numbers + // and rows. note that by the + // assertions above, we know that + // this particular shape function + // is primitive, so there is no + // question to which vector + // component the call of this + // function refers + return this->shape_values(shape_function_to_row_table[i], j); } + +template +inline +double +FEValuesBase::shape_value_component (const unsigned int i, + const unsigned int j, + const unsigned int component) const +{ + Assert (this->update_flags & update_values, + ExcAccessToUninitializedField()); + Assert (component < fe->n_components(), + ExcIndexRange(component, 0, fe->n_components())); + + // if this particulat shape + // function is primitive, then we + // can take a short-cut by checking + // whether the requested component + // is the only non-zero one (note + // that calling + // system_to_component_table only + // works if the shape function is + // primitive): + if (fe->is_primitive(i)) + { + if (component == fe->system_to_component_index(i).first) + return this->shape_values(i,j); + else + return 0; + } + else + { + // no, this shape function is + // not primitive. then we have + // to loop over its components + // and to find the + // corresponding row in the + // arrays of this + // object. before that check + // whether the shape function + // is non-zero at all within + // this component: + if (fe->get_nonzero_components(i)[component] == false) + return 0.; + + // count how many non-zero + // component the shape function + // has before the one we are + // looking for, and add this to + // the offset of the first + // non-zero component of this + // shape function in the arrays + // we index presently: + const unsigned int + row = (shape_function_to_row_table[i] + + + std::count (fe->get_nonzero_components(i).begin(), + fe->get_nonzero_components(i).begin()+component, + true)); + return this->shape_values(row, j); + }; +} + + + template inline const Tensor<1,dim> & FEValuesBase::shape_grad (const unsigned int i, const unsigned int j) const { - Assert (ishape_gradients.size(), ExcIndexRange (i, 0, shape_gradients.size())); Assert (jupdate_flags & update_gradients, + ExcAccessToUninitializedField()); + Assert (fe->is_primitive (i), + ExcShapeFunctionNotPrimitive(i)); + + // if the entire FE is primitive, + // then we can take a short-cut: + if (fe->is_primitive()) + return this->shape_gradients[i][j]; + else + // otherwise, use the mapping + // between shape function numbers + // and rows. note that by the + // assertions above, we know that + // this particular shape function + // is primitive, so there is no + // question to which vector + // component the call of this + // function refers + return this->shape_gradients[shape_function_to_row_table[i]][j]; }; +template +inline +Tensor<1,dim> +FEValuesBase::shape_grad_component (const unsigned int i, + const unsigned int j, + const unsigned int component) const +{ + Assert (this->update_flags & update_values, + ExcAccessToUninitializedField()); + Assert (component < fe->n_components(), + ExcIndexRange(component, 0, fe->n_components())); + + // if this particulat shape + // function is primitive, then we + // can take a short-cut by checking + // whether the requested component + // is the only non-zero one (note + // that calling + // system_to_component_table only + // works if the shape function is + // primitive): + if (fe->is_primitive(i)) + { + if (component == fe->system_to_component_index(i).first) + return this->shape_gradients[i][j]; + else + return Tensor<1,dim>(); + } + else + { + // no, this shape function is + // not primitive. then we have + // to loop over its components + // and to find the + // corresponding row in the + // arrays of this + // object. before that check + // whether the shape function + // is non-zero at all within + // this component: + if (fe->get_nonzero_components(i)[component] == false) + return Tensor<1,dim>(); + + // count how many non-zero + // component the shape function + // has before the one we are + // looking for, and add this to + // the offset of the first + // non-zero component of this + // shape function in the arrays + // we index presently: + const unsigned int + row = (shape_function_to_row_table[i] + + + std::count (fe->get_nonzero_components(i).begin(), + fe->get_nonzero_components(i).begin()+component, + true)); + return this->shape_gradients[row][j]; + }; +} + + + template inline const Tensor<2,dim> & FEValuesBase::shape_2nd_derivative (const unsigned int i, const unsigned int j) const { - Assert (ishape_2nd_derivatives.size(), ExcIndexRange (i, 0, shape_2nd_derivatives.size())); - Assert (jshape_2nd_derivatives[i].size(), ExcIndexRange (j, 0, shape_2nd_derivatives[i].size())); - Assert (update_flags & update_second_derivatives, ExcAccessToUninitializedField()); + Assert (this->update_flags & update_second_derivatives, + ExcAccessToUninitializedField()); + Assert (fe->is_primitive (i), + ExcShapeFunctionNotPrimitive(i)); + + // if the entire FE is primitive, + // then we can take a short-cut: + if (fe->is_primitive()) + return this->shape_2nd_derivatives[i][j]; + else + // otherwise, use the mapping + // between shape function numbers + // and rows. note that by the + // assertions above, we know that + // this particular shape function + // is primitive, so there is no + // question to which vector + // component the call of this + // function refers + return this->shape_2nd_derivatives[shape_function_to_row_table[i]][j]; +} + + - return shape_2nd_derivatives[i][j]; +template +inline +Tensor<2,dim> +FEValuesBase::shape_2nd_derivative_component (const unsigned int i, + const unsigned int j, + const unsigned int component) const +{ + Assert (this->update_flags & update_values, + ExcAccessToUninitializedField()); + Assert (component < fe->n_components(), + ExcIndexRange(component, 0, fe->n_components())); + + // if this particulat shape + // function is primitive, then we + // can take a short-cut by checking + // whether the requested component + // is the only non-zero one (note + // that calling + // system_to_component_table only + // works if the shape function is + // primitive): + if (fe->is_primitive(i)) + { + if (component == fe->system_to_component_index(i).first) + return this->shape_2nd_derivatives[i][j]; + else + return Tensor<2,dim>(); + } + else + { + // no, this shape function is + // not primitive. then we have + // to loop over its components + // and to find the + // corresponding row in the + // arrays of this + // object. before that check + // whether the shape function + // is non-zero at all within + // this component: + if (fe->get_nonzero_components(i)[component] == false) + return Tensor<2,dim>(); + + // count how many non-zero + // component the shape function + // has before the one we are + // looking for, and add this to + // the offset of the first + // non-zero component of this + // shape function in the arrays + // we index presently: + const unsigned int + row = (shape_function_to_row_table[i] + + + std::count (fe->get_nonzero_components(i).begin(), + fe->get_nonzero_components(i).begin()+component, + true)); + return this->shape_2nd_derivatives[row][j]; + }; } @@ -1065,11 +1505,11 @@ template const Point & FEFaceValuesBase::normal_vector (const unsigned int i) const { - Assert (inormal_vectors.size(), ExcIndexRange(i, 0, normal_vectors.size())); + Assert (this->update_flags & update_normal_vectors, FEValuesBase::ExcAccessToUninitializedField()); - return normal_vectors[i]; + return this->normal_vectors[i]; }; diff --git a/deal.II/deal.II/include/fe/mapping.h b/deal.II/deal.II/include/fe/mapping.h index 1240e8a725..41afb3822b 100644 --- a/deal.II/deal.II/include/fe/mapping.h +++ b/deal.II/deal.II/include/fe/mapping.h @@ -28,15 +28,16 @@ template class FEValues; template class FEFaceValues; template class FESubfaceValues; + /** - * Abstract basis class for mapping classes. + * Abstract base class for mapping classes. * * The interface for filling the tables of @ref{FEValues} is provided. * Everything else has to happen in derived classes. * * The following paragraph applies to the implementation of * @ref{FEValues}. Usage of the class is as follows: first, call the - * functionss @p{update_once} and @p{update_each} with the update + * functions @p{update_once} and @p{update_each} with the update * flags you need. This includes the flags needed by the * @ref{FiniteElement}. Then call @p{get_*_data} and with the or'd * results. This will initialize and return some internal data @@ -81,8 +82,41 @@ class Mapping : public Subscriptor const Point &p) const=0; /** - * Class for internal data of finite - * element and mapping objects. + * Base class for internal data + * of finite element and mapping + * objects. The internal + * mechanism is that upon + * construction of a @p{FEValues} + * objects, it asks the mapping + * and finite element classes + * that are to be used to + * allocate memory for their own + * purpose in which they may + * store data that only needs to + * be computed once. For example, + * most finite elements will + * store the values of the shape + * functions at the quadrature + * points in this object, since + * they do not change from cell + * to cell and only need to be + * computed once. Since different + * @p{FEValues} objects using + * different quadrature rules + * might access the same finite + * element object at the same + * time, it is necessary to + * create one such object per + * @p{FEValues} object. Ownership + * of this object is then + * transferred to the + * @p{FEValues} object, but a + * pointer to this object is + * passed to the finite element + * object every time it shall + * compute some data so that it + * has access to the precomputed + * values stored there. */ class InternalDataBase: public Subscriptor { diff --git a/deal.II/deal.II/source/dofs/dof_tools.cc b/deal.II/deal.II/source/dofs/dof_tools.cc index 85c9d652a8..c926fa205b 100644 --- a/deal.II/deal.II/source/dofs/dof_tools.cc +++ b/deal.II/deal.II/source/dofs/dof_tools.cc @@ -1715,7 +1715,7 @@ DoFTools::compute_intergrid_weights_1 (const DoFHandler &coars // the restriction of the fine // grid: const unsigned int coarse_dofs_per_cell_component - = coarse_fe.base_element(coarse_fe.component_to_base(coarse_component)).dofs_per_cell; + = coarse_fe.base_element(coarse_fe.component_to_base(coarse_component).first).dofs_per_cell; // Try to find out whether the @@ -1741,9 +1741,9 @@ DoFTools::compute_intergrid_weights_1 (const DoFHandler &coars ExcInvalidComponent (fine_component, fine_fe.n_components())); // check whether respective finite // elements are equal - Assert (coarse_fe.base_element (coarse_fe.component_to_base(coarse_component)) + Assert (coarse_fe.base_element (coarse_fe.component_to_base(coarse_component).first) == - fine_fe.base_element (fine_fe.component_to_base(fine_component)), + fine_fe.base_element (fine_fe.component_to_base(fine_component).first), ExcFiniteElementsDontMatch()); #ifdef DEBUG diff --git a/deal.II/deal.II/source/fe/fe.cc b/deal.II/deal.II/source/fe/fe.cc index 94bdcbcf25..50b0b86cf9 100644 --- a/deal.II/deal.II/source/fe/fe.cc +++ b/deal.II/deal.II/source/fe/fe.cc @@ -22,6 +22,10 @@ #include #include +#include +#include + + // if necessary try to work around a bug in the IBM xlC compiler #ifdef XLC_WORK_AROUND_STD_BUG using namespace std; @@ -105,19 +109,47 @@ FiniteElementBase::InternalDataBase::~InternalDataBase () template FiniteElementBase::FiniteElementBase (const FiniteElementData &fe_data, - const std::vector &restriction_is_additive_flags) + const std::vector &restriction_is_additive_flags, + const std::vector > &nonzero_components) : FiniteElementData (fe_data), system_to_component_table(dofs_per_cell), face_system_to_component_table(dofs_per_face), + system_to_base_table(dofs_per_cell), + face_system_to_base_table(dofs_per_face), component_to_system_table(components, std::vector(dofs_per_cell)), - face_component_to_system_table(components, std::vector(dofs_per_face)), - component_to_base_table (components, 0), - restriction_is_additive_flags(restriction_is_additive_flags) + face_component_to_system_table(components, std::vector(dofs_per_face)), + component_to_base_table (components, std::make_pair(0U, 0U)), + restriction_is_additive_flags(restriction_is_additive_flags), + nonzero_components (nonzero_components), + n_nonzero_components_table (compute_n_nonzero_components(nonzero_components)), + cached_primitivity (std::find_if (n_nonzero_components_table.begin(), + n_nonzero_components_table.end(), + std::bind2nd(std::not_equal_to(), + 1U)) + == + n_nonzero_components_table.end()) { - Assert(restriction_is_additive_flags.size()==fe_data.components, - ExcDimensionMismatch(restriction_is_additive_flags.size(),fe_data.components)); - + Assert (restriction_is_additive_flags.size()==fe_data.components, + ExcDimensionMismatch(restriction_is_additive_flags.size(), + fe_data.components)); + Assert (nonzero_components.size() == dofs_per_cell, + ExcInternalError()); + for (unsigned int i=0; i= 1, + ExcInternalError()); + Assert (n_nonzero_components_table[i] >= 1, + ExcInternalError()); + Assert (n_nonzero_components_table[i] <= n_components(), + ExcInternalError()); + }; + for (unsigned int i=0; i::children_per_cell; ++i) { restriction[i].reinit (dofs_per_cell, dofs_per_cell); @@ -159,11 +191,13 @@ FiniteElementBase::FiniteElementBase (const FiniteElementData &fe_data for (unsigned int j=0 ; j(0,j); + system_to_base_table[j] = std::make_pair(std::make_pair(0U,0U),j); component_to_system_table[0][j] = j; } for (unsigned int j=0 ; j(0,j); + face_system_to_base_table[j] = std::make_pair(std::make_pair(0U,0U),j); face_component_to_system_table[0][j] = j; } }; @@ -180,6 +214,18 @@ FiniteElementBase::shape_value (const unsigned int, +template +double +FiniteElementBase::shape_value_component (const unsigned int, + const Point &, + const unsigned int) const +{ + AssertThrow(false, ExcUnitShapeValuesDoNotExist()); + return 0.; +} + + + template Tensor<1,dim> FiniteElementBase::shape_grad (const unsigned int, @@ -191,6 +237,18 @@ FiniteElementBase::shape_grad (const unsigned int, +template +Tensor<1,dim> +FiniteElementBase::shape_grad_component (const unsigned int, + const Point &, + const unsigned int) const +{ + AssertThrow(false, ExcUnitShapeValuesDoNotExist()); + return Tensor<1,dim> (); +} + + + template Tensor<2,dim> FiniteElementBase::shape_grad_grad (const unsigned int, @@ -201,6 +259,18 @@ FiniteElementBase::shape_grad_grad (const unsigned int, } + +template +Tensor<2,dim> +FiniteElementBase::shape_grad_grad_component (const unsigned int, + const Point &, + const unsigned int) const +{ + AssertThrow(false, ExcUnitShapeValuesDoNotExist()); + return Tensor<2,dim> (); +} + + template const FullMatrix & FiniteElementBase::restrict (const unsigned int child) const @@ -316,10 +386,14 @@ FiniteElementBase::memory_consumption () const MemoryConsumption::memory_consumption (interface_constraints) + MemoryConsumption::memory_consumption (system_to_component_table) + MemoryConsumption::memory_consumption (face_system_to_component_table) + + MemoryConsumption::memory_consumption (system_to_base_table) + + MemoryConsumption::memory_consumption (face_system_to_base_table) + MemoryConsumption::memory_consumption (component_to_system_table) + MemoryConsumption::memory_consumption (face_component_to_system_table) + MemoryConsumption::memory_consumption (component_to_base_table) + - MemoryConsumption::memory_consumption (restriction_is_additive_flags)); + MemoryConsumption::memory_consumption (restriction_is_additive_flags) + + MemoryConsumption::memory_consumption (nonzero_components) + + MemoryConsumption::memory_consumption (n_nonzero_components_table)); }; @@ -402,13 +476,31 @@ compute_2nd (const Mapping &mapping, } + +template +std::vector +FiniteElementBase:: +compute_n_nonzero_components (const std::vector > &nonzero_components) +{ + std::vector retval (nonzero_components.size()); + for (unsigned int i=0; i FiniteElement::FiniteElement (const FiniteElementData &fe_data, - const std::vector &restriction_is_additive_flags) : + const std::vector &restriction_is_additive_flags, + const std::vector > &nonzero_components) : FiniteElementBase (fe_data, - restriction_is_additive_flags) + restriction_is_additive_flags, + nonzero_components) {} diff --git a/deal.II/deal.II/source/fe/fe_dgp.cc b/deal.II/deal.II/source/fe/fe_dgp.cc index 2e6db63f43..12678e7777 100644 --- a/deal.II/deal.II/source/fe/fe_dgp.cc +++ b/deal.II/deal.II/source/fe/fe_dgp.cc @@ -12,8 +12,6 @@ //---------------------------------------------------------------- #include -#include -#include #include #include #include @@ -28,16 +26,12 @@ template FE_DGP::FE_DGP (unsigned int degree) : FiniteElement (FiniteElementData(get_dpo_vector(degree),1), - std::vector(1,true)), + std::vector(1,true), + std::vector >(FiniteElementData(get_dpo_vector(degree),1).dofs_per_cell, + std::vector(1,true))), degree(degree), - poly(0) + polynomial_space (Legendre::generate_complete_basis(degree)) { - // create array of Legendre polynomials - std::vector > v; - for (unsigned int i=0;i<=degree;++i) - v.push_back(Legendre(i)); - poly = new PolynomialSpace (v); - // if defined, copy over matrices // from precomputed arrays if ((degree < Matrices::n_embedding_matrices) && @@ -47,6 +41,7 @@ FE_DGP::FE_DGP (unsigned int degree) else for (unsigned int i=0; i::children_per_cell;++i) prolongation[i].reinit(0,0); + // // same as above: copy over matrix // // from predefined values and // // generate all others by rotation @@ -64,31 +59,34 @@ FE_DGP::FE_DGP (unsigned int degree) template -FE_DGP::~FE_DGP () +FiniteElement * +FE_DGP::clone() const { - // delete poly member and set it to - // zero to prevent accidental use - delete poly; - poly = 0; + return new FE_DGP(degree); } template -FiniteElement * -FE_DGP::clone() const +double +FE_DGP::shape_value (const unsigned int i, + const Point &p) const { - return new FE_DGP(degree); + Assert (i double -FE_DGP::shape_value (const unsigned int i, - const Point &p) const +FE_DGP::shape_value_component (const unsigned int i, + const Point &p, + const unsigned int component) const { - return poly->compute_value(i, p); + Assert (i FE_DGP::shape_grad (const unsigned int i, const Point &p) const { - return poly->compute_grad(i, p); + Assert (i +Tensor<1,dim> +FE_DGP::shape_grad_component (const unsigned int i, + const Point &p, + const unsigned int component) const +{ + Assert (i FE_DGP::shape_grad_grad (const unsigned int i, const Point &p) const { - return poly->compute_grad_grad(i, p); + Assert (i +Tensor<2,dim> +FE_DGP::shape_grad_grad_component (const unsigned int i, + const Point &p, + const unsigned int component) const +{ + Assert (i UpdateFlags FE_DGP::update_once (const UpdateFlags flags) const { - UpdateFlags out = update_default; - - if (flags & update_values) - out |= update_values; - - return out; + // for this kind of elements, only + // the values can be precomputed + // once and for all. set this flag + // if the values are requested at + // all + return (update_default | (flags & update_values)); } @@ -171,11 +196,8 @@ FE_DGP::get_data (const UpdateFlags update_flags, const Mapping &mapping, const Quadrature &quadrature) const { + // generate a new data object InternalData* data = new InternalData; - std::vector values(0); - std::vector > grads(0); - std::vector > grad_grads(0); - // check what needs to be // initialized only once and what // on every cell/face/subface we @@ -186,7 +208,15 @@ FE_DGP::get_data (const UpdateFlags update_flags, const UpdateFlags flags(data->update_flags); const unsigned int n_q_points = quadrature.n_quadrature_points; + + // have some scratch arrays + std::vector values(0); + std::vector > grads(0); + std::vector > grad_grads(0); + // initialize fields only if really + // necessary. otherwise, don't + // allocate memory if (flags & update_values) { values.resize (dofs_per_cell); @@ -208,11 +238,18 @@ FE_DGP::get_data (const UpdateFlags update_flags, if (flags & update_second_derivatives) data->initialize_2nd (this, mapping, quadrature); - + // next already fill those fields + // of which we have information by + // now. note that the shape + // gradients are only those on the + // unit cell, and need to be + // transformed when visiting an + // actual cell if (flags & (update_values | update_gradients)) for (unsigned int i=0; icompute(quadrature.point(i), values, grads, grad_grads); + polynomial_space.compute(quadrature.point(i), + values, grads, grad_grads); for (unsigned int k=0; k::fill_fe_values (const Mapping &mapping, for (unsigned int k=0; k FE_DGQ::FE_DGQ (unsigned int degree) : FiniteElement (FiniteElementData(get_dpo_vector(degree),1), - std::vector(1,true)), + std::vector(1,true), + std::vector >(FiniteElementData(get_dpo_vector(degree),1).dofs_per_cell, + std::vector(1,true))), degree(degree), - poly(0) + polynomial_space (LagrangeEquidistant::generate_complete_basis(degree)) { - if (degree==0) - { - // create constant polynomial - std::vector > - v(1, Polynomial (std::vector (1,1.))); - poly = new TensorProductPolynomials (v); - } - else - { - // create array of Lagrange polynomials - std::vector v; - for (unsigned int i=0;i<=degree;++i) - v.push_back(LagrangeEquidistant(degree,i)); - poly = new TensorProductPolynomials (v); - } - // generate permutation/rotation // index sets to generate some // matrices from others @@ -194,31 +180,34 @@ FE_DGQ::FE_DGQ (unsigned int degree) template -FE_DGQ::~FE_DGQ () +FiniteElement * +FE_DGQ::clone() const { - // delete poly member and set it to - // zero to prevent accidental use - delete poly; - poly = 0; + return new FE_DGQ(degree); } template -FiniteElement * -FE_DGQ::clone() const +double +FE_DGQ::shape_value (const unsigned int i, + const Point &p) const { - return new FE_DGQ(degree); + Assert (i double -FE_DGQ::shape_value (const unsigned int i, - const Point &p) const +FE_DGQ::shape_value_component (const unsigned int i, + const Point &p, + const unsigned int component) const { - return poly->compute_value(i, p); + Assert (i FE_DGQ::shape_grad (const unsigned int i, const Point &p) const { - return poly->compute_grad(i, p); + Assert (i +Tensor<1,dim> +FE_DGQ::shape_grad_component (const unsigned int i, + const Point &p, + const unsigned int component) const +{ + Assert (i FE_DGQ::shape_grad_grad (const unsigned int i, const Point &p) const { - return poly->compute_grad_grad(i, p); + Assert (i +Tensor<2,dim> +FE_DGQ::shape_grad_grad_component (const unsigned int i, + const Point &p, + const unsigned int component) const +{ + Assert (i UpdateFlags FE_DGQ::update_once (const UpdateFlags flags) const { - UpdateFlags out = update_default; - - if (flags & update_values) - out |= update_values; - - return out; + // for this kind of elements, only + // the values can be precomputed + // once and for all. set this flag + // if the values are requested at + // all + return (update_default | (flags & update_values)); } @@ -379,10 +395,8 @@ FE_DGQ::get_data (const UpdateFlags update_flags, const Mapping &mapping, const Quadrature &quadrature) const { + // generate a new data object InternalData* data = new InternalData; - std::vector values(0); - std::vector > grads(0); - std::vector > grad_grads(0); // check what needs to be // initialized only once and what @@ -394,7 +408,15 @@ FE_DGQ::get_data (const UpdateFlags update_flags, const UpdateFlags flags(data->update_flags); const unsigned int n_q_points = quadrature.n_quadrature_points; - + + // have some scratch arrays + std::vector values(0); + std::vector > grads(0); + std::vector > grad_grads(0); + + // initialize fields only if really + // necessary. otherwise, don't + // allocate memory if (flags & update_values) { values.resize (dofs_per_cell); @@ -416,11 +438,18 @@ FE_DGQ::get_data (const UpdateFlags update_flags, if (flags & update_second_derivatives) data->initialize_2nd (this, mapping, quadrature); - + // next already fill those fields + // of which we have information by + // now. note that the shape + // gradients are only those on the + // unit cell, and need to be + // transformed when visiting an + // actual cell if (flags & (update_values | update_gradients)) for (unsigned int i=0; icompute(quadrature.point(i), values, grads, grad_grads); + polynomial_space.compute(quadrature.point(i), + values, grads, grad_grads); for (unsigned int k=0; k::fill_fe_values (const Mapping &mapping, for (unsigned int k=0; k::has_support_on_face (const unsigned int shape_index, const unsigned int face_index) const { + Assert (shape_index < dofs_per_cell, + ExcIndexRange (shape_index, 0, dofs_per_cell)); + Assert (face_index < GeometryInfo::faces_per_cell, + ExcIndexRange (face_index, 0, GeometryInfo::faces_per_cell)); unsigned int n = degree+1; unsigned int n2 = n*n; diff --git a/deal.II/deal.II/source/fe/fe_nedelec.cc b/deal.II/deal.II/source/fe/fe_nedelec.cc new file mode 100644 index 0000000000..2552e4824b --- /dev/null +++ b/deal.II/deal.II/source/fe/fe_nedelec.cc @@ -0,0 +1,863 @@ +//---------------------------------------------------------------- +// $Id$ +// Version: $Name$ +// +// Copyright (C) 2002 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 + + +template +FE_Nedelec::FE_Nedelec (const unsigned int degree) + : + FiniteElement (FiniteElementData(get_dpo_vector(degree), + dim), +//TODO: I'd think this element is actually additive in the restriction + std::vector (dim,false), + std::vector >(FiniteElementData(get_dpo_vector(degree),dim).dofs_per_cell, + std::vector(dim,true))), + degree(degree) +{ + Assert (dim >= 2, ExcNotUsefulInThisDimension()); + + // copy constraint matrices if they + // are defined. otherwise set them + // to invalid size + if (degree::children_per_cell; ++c) + prolongation[c].fill (Matrices::embedding[degree-1][c]); + else + for (unsigned int i=0; i::children_per_cell;++i) + prolongation[i].reinit(0,0); + + // then fill restriction + // matrices. they are hardcoded for + // the first few elements + switch (dim) + { + case 2: // 2d + { + switch (degree) + { + case 1: + { + // DoF on bottom line + // of coarse cell will + // be mean value of + // bottom DoFs on the + // two adjacent child + // cells + restriction[0](0,0) = 0.5; + restriction[1](0,0) = 0.5; + // same for other DoFs + restriction[1](1,1) = 0.5; + restriction[2](1,1) = 0.5; + + restriction[2](2,2) = 0.5; + restriction[3](2,2) = 0.5; + + restriction[3](3,3) = 0.5; + restriction[0](3,3) = 0.5; + + break; + }; + + default: + { + // in case we don't + // have the matrices + // (yet), set them to + // impossible + // values. this does + // not prevent the use + // of this FE, but will + // prevent the use of + // these matrices + for (unsigned int i=0; + i::children_per_cell; + ++i) + restriction[i].reinit(0,0); + }; + }; + + break; + }; + + + case 3: // 3d + { + switch (degree) + { + case 1: + { + // same principle as in + // 2d + restriction[0](0,0) = 0.5; + restriction[1](0,0) = 0.5; + + restriction[1](1,1) = 0.5; + restriction[2](1,1) = 0.5; + + restriction[2](2,2) = 0.5; + restriction[3](2,2) = 0.5; + + restriction[3](3,3) = 0.5; + restriction[0](3,3) = 0.5; + + restriction[4](4,4) = 0.5; + restriction[5](4,4) = 0.5; + + restriction[5](5,5) = 0.5; + restriction[6](5,5) = 0.5; + + restriction[6](6,6) = 0.5; + restriction[7](6,6) = 0.5; + + restriction[7](7,7) = 0.5; + restriction[4](7,7) = 0.5; + + + restriction[1](8,8) = 0.5; + restriction[5](8,8) = 0.5; + + restriction[2](9,9) = 0.5; + restriction[6](9,9) = 0.5; + + restriction[3](10,10) = 0.5; + restriction[7](10,10) = 0.5; + + restriction[0](11,11) = 0.5; + restriction[5](11,11) = 0.5; + + break; + }; + + default: + { + // in case we don't + // have the matrices + // (yet), set them to + // impossible + // values. this does + // not prevent the use + // of this FE, but will + // prevent the use of + // these matrices + for (unsigned int i=0; + i::children_per_cell; + ++i) + restriction[i].reinit(0,0); + }; + }; + + break; + }; + + default: + Assert (false,ExcNotImplemented()); + } + + // finally fill in support points + // on cell and face + initialize_unit_support_points (); + initialize_unit_face_support_points (); +}; + + + +template +FiniteElement * +FE_Nedelec::clone() const +{ + return new FE_Nedelec(degree); +} + + + +template +double +FE_Nedelec::shape_value_component (const unsigned int i, + const Point &p, + const unsigned int component) const +{ + Assert (i=8) && (component != 1))) + return 0; + + // now we know that the + // only non-zero + // component is + // requested: +//TODO[Anna]: check + const double x = p(0), + y = p(1), + z = p(2); + switch (i) + { + case 0: return (1-y)*(1-z); + case 2: return (1-y)*z; + case 1: return x*(1-y); + case 3: return (1-x)*(1-y); + + case 4: return y*(1-z); + case 6: return y*z; + case 5: return x*y; + case 7: return (1-x)*y; + + case 8: return (1-x)*(1-z); + case 9: return x*(1-z); + case 10: return x*z; + case 11: return (1-x)*z; + default: + Assert (false, ExcInternalError()); + return 0; + }; + }; + + // no other degrees + // implemented + default: + Assert (false, ExcNotImplemented()); + }; + }; + + // presently no other space + // dimension implemented + default: + Assert (false, ExcNotImplemented()); + }; + + return 0; +} + + + +template +Tensor<1,dim> +FE_Nedelec::shape_grad_component (const unsigned int i, + const Point &/*p*/, + const unsigned int component) const +{ + Assert (i(unit_gradients[i][component]); + }; + + // no other degrees + // implemented + default: + Assert (false, ExcNotImplemented()); + }; + }; + + // presently no other space + // dimension implemented + default: + Assert (false, ExcNotImplemented()); + }; + + return Tensor<1,dim>(); +} + + + +template +Tensor<2,dim> +FE_Nedelec::shape_grad_grad_component (const unsigned int i, + const Point &/*p*/, + const unsigned int component) const +{ + Assert (i(); + }; + + // no other degrees + // implemented + default: + Assert (false, ExcNotImplemented()); + }; + }; + + case 3: // 3D + { + switch (degree) + { + // first order Nedelec + // elements. their second + // derivatives on the + // unit cell are zero + case 1: + { + return Tensor<2,dim>(); + }; + + // no other degrees + // implemented + default: + Assert (false, ExcNotImplemented()); + }; + }; + + + // presently no other space + // dimension implemented + default: + Assert (false, ExcNotImplemented()); + }; + + return Tensor<2,dim>(); +} + + +//---------------------------------------------------------------------- +// Auxiliary functions +//---------------------------------------------------------------------- + + + +template +void FE_Nedelec::initialize_unit_support_points () +{ +//TODO: fix for higher orders. correct now for lowest order, all dimensions +// is this correct? all DoFs on lines, none on faces or bubbles? + + // all degrees of freedom are on + // edges, and their order is the + // same as the edges themselves + unit_support_points.resize(GeometryInfo::lines_per_cell * degree); + unsigned int index = 0; + for (unsigned int line=0; line::lines_per_cell; ++line) + { + const unsigned int + vertex_index_0 = GeometryInfo::vertices_adjacent_to_line(line,0), + vertex_index_1 = GeometryInfo::vertices_adjacent_to_line(line,1); + + const Point + vertex_0 = GeometryInfo::unit_cell_vertex(vertex_index_0), + vertex_1 = GeometryInfo::unit_cell_vertex(vertex_index_1); + + // place dofs equispaced + // between the vertices of each + // line + for (unsigned int d=0; d +void FE_Nedelec<1>::initialize_unit_face_support_points () +{ + // no faces in 1d, so nothing to do +}; + +#endif + + +template +void FE_Nedelec::initialize_unit_face_support_points () +{ +//TODO: fix for higher orders. correct now for lowest order, all dimensions +// is this correct? all DoFs on lines, none on faces or bubbles? + // do this the same as above, but + // for one dimension less + unit_face_support_points.resize(GeometryInfo::lines_per_cell * degree); + unsigned int index = 0; + for (unsigned int line=0; line::lines_per_cell; ++line) + { + const unsigned int + vertex_index_0 = GeometryInfo::vertices_adjacent_to_line(line,0), + vertex_index_1 = GeometryInfo::vertices_adjacent_to_line(line,1); + + const Point + vertex_0 = GeometryInfo::unit_cell_vertex(vertex_index_0), + vertex_1 = GeometryInfo::unit_cell_vertex(vertex_index_1); + + // place dofs equispaced + // between the vertices of each + // line + for (unsigned int d=0; d +std::vector +FE_Nedelec::get_dpo_vector(const unsigned int degree) +{ +//TODO: fix for higher orders. correct now for lowest order, all dimensions + std::vector dpo(dim+1, 0); +// can this be done in a dimension independent and degree independent way? +// if DoFs are located only on lines, the the following is the correct way + + // put all degrees of freedom on + // the lines, and in particular + // @p{degree} DoFs per line: + dpo[1] = degree; + return dpo; +} + + + +template +UpdateFlags +FE_Nedelec::update_once (const UpdateFlags flags) const +{ +//TODO: think about what this actually means here??? + // for this kind of elements, only + // the values can be precomputed + // once and for all. set this flag + // if the values are requested at + // all + return (update_default | (flags & update_values)); +} + + + +template +UpdateFlags +FE_Nedelec::update_each (const UpdateFlags flags) const +{ +//TODO: think about what this actually means here??? + + UpdateFlags out = update_default; + + if (flags & update_gradients) + out |= update_gradients | update_covariant_transformation; + if (flags & update_second_derivatives) + out |= update_second_derivatives | update_covariant_transformation; + + return out; +} + + + +//---------------------------------------------------------------------- +// Data field initialization +//---------------------------------------------------------------------- + +template +typename Mapping::InternalDataBase * +FE_Nedelec::get_data (const UpdateFlags /*update_flags*/, + const Mapping &/*mapping*/, + const Quadrature &/*quadrature*/) const +{ + return 0; +//TODO +// // generate a new data object and +// // initialize some fields +// InternalData* data = new InternalData; + +// // check what needs to be +// // initialized only once and what +// // on every cell/face/subface we +// // visit +// 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.n_quadrature_points; + +// // some scratch arrays +// std::vector values(0); +// std::vector > grads(0); +// std::vector > grad_grads(0); + +// // initialize fields only if really +// // necessary. otherwise, don't +// // allocate memory +// if (flags & update_values) +// { +// values.resize (dofs_per_cell); +// data->shape_values.resize(dofs_per_cell, +// std::vector(n_q_points)); +// } + +// if (flags & update_gradients) +// { +// grads.resize (dofs_per_cell); +// data->shape_gradients.resize(dofs_per_cell, +// std::vector >(n_q_points)); +// } + +// // if second derivatives through +// // finite differencing is required, +// // then initialize some objects for +// // that +// if (flags & update_second_derivatives) +// data->initialize_2nd (this, mapping, quadrature); + +// // next already fill those fields +// // of which we have information by +// // now. note that the shape +// // gradients are only those on the +// // unit cell, and need to be +// // transformed when visiting an +// // actual cell +// if (flags & (update_values | update_gradients)) +// for (unsigned int i=0; ishape_values[renumber[k]][i] = values[k]; + +// if (flags & update_gradients) +// for (unsigned int k=0; kshape_gradients[renumber[k]][i] = grads[k]; +// } +// return data; +} + + + + +//---------------------------------------------------------------------- +// Fill data of FEValues +//---------------------------------------------------------------------- + +template +void +FE_Nedelec::fill_fe_values (const Mapping &/*mapping*/, + const typename DoFHandler::cell_iterator &/*cell*/, + const Quadrature &/*quadrature*/, + typename Mapping::InternalDataBase &/*mapping_data*/, + typename Mapping::InternalDataBase &/*fedata*/, + FEValuesData &/*data*/) const +{ +//TODO!! +// // convert data object to internal +// // data for this class. fails with +// // an exception if that is not +// // possible +// InternalData &fe_data = dynamic_cast (fedata); + +// const UpdateFlags flags(fe_data.current_update_flags()); + +// for (unsigned int k=0; k +void +FE_Nedelec::fill_fe_face_values (const Mapping &/*mapping*/, + const typename DoFHandler::cell_iterator &/*cell*/, + const unsigned int /*face*/, + const Quadrature &/*quadrature*/, + typename Mapping::InternalDataBase &/*mapping_data*/, + typename Mapping::InternalDataBase &/*fedata*/, + FEValuesData &/*data*/) const +{ +//TODO!! +// // convert data object to internal +// // data for this class. fails with +// // an exception if that is not +// // possible +// InternalData &fe_data = dynamic_cast (fedata); + +// // offset determines which data set +// // to take (all data sets for all +// // faces are stored contiguously) +// const unsigned int offset = face * quadrature.n_quadrature_points; + +// const UpdateFlags flags(fe_data.update_once | fe_data.update_each); + +// for (unsigned int k=0; k +void +FE_Nedelec::fill_fe_subface_values (const Mapping &/*mapping*/, + const typename DoFHandler::cell_iterator &/*cell*/, + const unsigned int /*face*/, + const unsigned int /*subface*/, + const Quadrature &/*quadrature*/, + typename Mapping::InternalDataBase &/*mapping_data*/, + typename Mapping::InternalDataBase &/*fedata*/, + FEValuesData &/*data*/) const +{ +//TODO!! +// // convert data object to internal +// // data for this class. fails with +// // an exception if that is not +// // possible +// InternalData &fe_data = dynamic_cast (fedata); + +// // offset determines which data set +// // to take (all data sets for all +// // sub-faces are stored contiguously) +// const unsigned int offset = (face * GeometryInfo::subfaces_per_face + subface) +// * quadrature.n_quadrature_points; + +// const UpdateFlags flags(fe_data.update_once | fe_data.update_each); + +// for (unsigned int k=0; k +unsigned int +FE_Nedelec::n_base_elements () const +{ + return 1; +}; + + + +template +const FiniteElement & +FE_Nedelec::base_element (const unsigned int index) const +{ + Assert (index==0, ExcIndexRange(index, 0, 1)); + return *this; +}; + + + +template +bool +FE_Nedelec::has_support_on_face (const unsigned int shape_index, + const unsigned int face_index) const +{ + Assert (shape_index < dofs_per_cell, + ExcIndexRange (shape_index, 0, dofs_per_cell)); + Assert (face_index < GeometryInfo::faces_per_cell, + ExcIndexRange (face_index, 0, GeometryInfo::faces_per_cell)); + +//TODO: fix for higher orders. correct now for lowest order, all dimensions +//TODO!! +// can this be done in a way that is dimension and degree independent? + + // all degrees of freedom are on + // lines, so also on a face. the + // question is whether it has + // support on this particular face + Assert (false, ExcNotImplemented()); + return true; +} + + + +template +unsigned int +FE_Nedelec::memory_consumption () const +{ + Assert (false, ExcNotImplemented ()); + return 0; +} + + + +template +unsigned int +FE_Nedelec::get_degree () const +{ + return degree; +}; + + + +template class FE_Nedelec; diff --git a/deal.II/deal.II/source/fe/fe_nedelec_1d.cc b/deal.II/deal.II/source/fe/fe_nedelec_1d.cc new file mode 100644 index 0000000000..10cf287383 --- /dev/null +++ b/deal.II/deal.II/source/fe/fe_nedelec_1d.cc @@ -0,0 +1,53 @@ +//---------------------------------------------------------------- +// $Id$ +// Version: $Name$ +// +// Copyright (C) 2002 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. +// +//---------------------------------------------------------------- + + + +// only compile this file if in 1d. note that Nedelec elemets do not +// make much sense in 1d, so this file only contains dummy +// implementations to avoid linker errors due to missing symbols +#if deal_II_dimension == 1 + + +#include + + +template <> +const double * const +FE_Nedelec<1>::Matrices::embedding[][GeometryInfo<1>::children_per_cell] = +{}; + + +template <> +const unsigned int +FE_Nedelec<1>::Matrices::n_embedding_matrices = 0; + + + +// No constraints in 1d +template <> +const unsigned int +FE_Nedelec<1>::Matrices::n_constraint_matrices = 0; + + +template <> +const double * const +FE_Nedelec<1>::Matrices::constraint_matrices[] = {}; + + +#else // #if deal_II_dimension +// On gcc2.95 on Alpha OSF1, the native assembler does not like empty +// files, so provide some dummy code +namespace { void dummy () {}; }; +#endif // #if deal_II_dimension == 1 + diff --git a/deal.II/deal.II/source/fe/fe_nedelec_2d.cc b/deal.II/deal.II/source/fe/fe_nedelec_2d.cc new file mode 100644 index 0000000000..fd7f8a5775 --- /dev/null +++ b/deal.II/deal.II/source/fe/fe_nedelec_2d.cc @@ -0,0 +1,118 @@ +//---------------------------------------------------------------- +// $Id$ +// Version: $Name$ +// +// Copyright (C) 2001, 2002 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. +// +//---------------------------------------------------------------- + + +// only compile this file if in 2d +#if deal_II_dimension == 2 + + +#include + +// Transfer matrices for finite elements: have one matrix for each of +// the four child cells which tells us how the degrees of freedom on +// the child cell are obtained from the degrees of freedom on the +// mother cell +namespace FE_Nedelec_2d +{ + static const double q1_into_q1_refined_0[] = + { + 1., 0, 0, 0, + 0, 0.5,0, 0.5, + 0.5, 0, 0.5,0, + 0, 0, 0, 1 + }; + + static const double q1_into_q1_refined_1[] = + { +//TODO[Anna] +// 13.5/27., 13.5/27., 0., 0., +// 0., 1., 0., 0., +// 0., 13.5/27., 13.5/27., 0., +// 6.75/27., 6.75/27., 6.75/27., 6.75/27., + }; + + static const double q1_into_q1_refined_2[] = + { +// 6.75/27., 6.75/27., 6.75/27., 6.75/27., +// 0., 13.5/27., 13.5/27., 0., +// 0., 0., 1., 0., +// 0., 0., 13.5/27., 13.5/27., + }; + + static const double q1_into_q1_refined_3[] = + { +// 13.5/27., 0., 0., 13.5/27., +// 6.75/27., 6.75/27., 6.75/27., 6.75/27., +// 0., 0., 13.5/27., 13.5/27., +// 0., 0., 0., 1., + }; +}; // namespace FE_Nedelec_2d + + +// embedding matrices + +template <> +const double * const +FE_Nedelec<2>::Matrices::embedding[][GeometryInfo<2>::children_per_cell] = +{ + { FE_Nedelec_2d::q1_into_q1_refined_0, FE_Nedelec_2d::q1_into_q1_refined_1, + FE_Nedelec_2d::q1_into_q1_refined_2, FE_Nedelec_2d::q1_into_q1_refined_3 } +}; + + +template <> +const unsigned int +FE_Nedelec<2>::Matrices::n_embedding_matrices + = sizeof(FE_Nedelec<2>::Matrices::embedding) / + sizeof(FE_Nedelec<2>::Matrices::embedding[0]); + + +// Constraint matrices: how do the new value on child faces depend on +// the values on the mother face if that face has a hanging node +namespace FE_Nedelec_2d +{ + static const double constraint_q1[] = + { + // the function is constant + // along each edge, so each + // degree of freedom on the + // refined edge has the same + // value as that on the + // coarse edge + 1., 1. + }; + +}; + + +template <> +const double * const +FE_Nedelec<2>::Matrices::constraint_matrices[] = +{ + FE_Nedelec_2d::constraint_q1 +}; + + +template <> +const unsigned int +FE_Nedelec<2>::Matrices::n_constraint_matrices + = sizeof(FE_Nedelec<2>::Matrices::constraint_matrices) / + sizeof(FE_Nedelec<2>::Matrices::constraint_matrices[0]); + + + +#else // #if deal_II_dimension +// On gcc2.95 on Alpha OSF1, the native assembler does not like empty +// files, so provide some dummy code +namespace { void dummy () {}; }; +#endif // #if deal_II_dimension == 2 diff --git a/deal.II/deal.II/source/fe/fe_nedelec_3d.cc b/deal.II/deal.II/source/fe/fe_nedelec_3d.cc new file mode 100644 index 0000000000..425fde8525 --- /dev/null +++ b/deal.II/deal.II/source/fe/fe_nedelec_3d.cc @@ -0,0 +1,183 @@ +//---------------------------------------------------------------- +// $Id$ +// Version: $Name$ +// +// Copyright (C) 2001, 2002 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. +// +//---------------------------------------------------------------- + +// Transfer matrices for finite elements + + +// only compile this file if in 3d +#if deal_II_dimension == 3 + +#include + +namespace FE_Nedelec_3d +{ + static const double q1_into_q1_refined_0[] = + { +// 1., 0., 0., 0., 0., 0., 0., 0., +// 13.5/27., 13.5/27., 0., 0., 0., 0., 0., 0., +// 6.75/27., 6.75/27., 6.75/27., 6.75/27., 0., 0., 0., 0., +// 13.5/27., 0., 0., 13.5/27., 0., 0., 0., 0., +// 13.5/27., 0., 0., 0., 13.5/27., 0., 0., 0., +// 6.75/27., 6.75/27., 0., 0., 6.75/27., 6.75/27., 0., 0., +// 3.375/27., 3.375/27., 3.375/27., 3.375/27., 3.375/27., 3.375/27., 3.375/27., 3.375/27., +// 6.75/27., 0., 0., 6.75/27., 6.75/27., 0., 0., 6.75/27., + }; + + static const double q1_into_q1_refined_1[] = + { +// 13.5/27., 13.5/27., 0., 0., 0., 0., 0., 0., +// 0., 1., 0., 0., 0., 0., 0., 0., +// 0., 13.5/27., 13.5/27., 0., 0., 0., 0., 0., +// 6.75/27., 6.75/27., 6.75/27., 6.75/27., 0., 0., 0., 0., +// 6.75/27., 6.75/27., 0., 0., 6.75/27., 6.75/27., 0., 0., +// 0., 13.5/27., 0., 0., 0., 13.5/27., 0., 0., +// 0., 6.75/27., 6.75/27., 0., 0., 6.75/27., 6.75/27., 0., +// 3.375/27., 3.375/27., 3.375/27., 3.375/27., 3.375/27., 3.375/27., 3.375/27., 3.375/27., + }; + + static const double q1_into_q1_refined_2[] = + { +// 6.75/27., 6.75/27., 6.75/27., 6.75/27., 0., 0., 0., 0., +// 0., 13.5/27., 13.5/27., 0., 0., 0., 0., 0., +// 0., 0., 1., 0., 0., 0., 0., 0., +// 0., 0., 13.5/27., 13.5/27., 0., 0., 0., 0., +// 3.375/27., 3.375/27., 3.375/27., 3.375/27., 3.375/27., 3.375/27., 3.375/27., 3.375/27., +// 0., 6.75/27., 6.75/27., 0., 0., 6.75/27., 6.75/27., 0., +// 0., 0., 13.5/27., 0., 0., 0., 13.5/27., 0., +// 0., 0., 6.75/27., 6.75/27., 0., 0., 6.75/27., 6.75/27., + }; + + static const double q1_into_q1_refined_3[] = + { +// 13.5/27., 0., 0., 13.5/27., 0., 0., 0., 0., +// 6.75/27., 6.75/27., 6.75/27., 6.75/27., 0., 0., 0., 0., +// 0., 0., 13.5/27., 13.5/27., 0., 0., 0., 0., +// 0., 0., 0., 1., 0., 0., 0., 0., +// 6.75/27., 0., 0., 6.75/27., 6.75/27., 0., 0., 6.75/27., +// 3.375/27., 3.375/27., 3.375/27., 3.375/27., 3.375/27., 3.375/27., 3.375/27., 3.375/27., +// 0., 0., 6.75/27., 6.75/27., 0., 0., 6.75/27., 6.75/27., +// 0., 0., 0., 13.5/27., 0., 0., 0., 13.5/27., + }; + + static const double q1_into_q1_refined_4[] = + { +// 13.5/27., 0., 0., 0., 13.5/27., 0., 0., 0., +// 6.75/27., 6.75/27., 0., 0., 6.75/27., 6.75/27., 0., 0., +// 3.375/27., 3.375/27., 3.375/27., 3.375/27., 3.375/27., 3.375/27., 3.375/27., 3.375/27., +// 6.75/27., 0., 0., 6.75/27., 6.75/27., 0., 0., 6.75/27., +// 0., 0., 0., 0., 1., 0., 0., 0., +// 0., 0., 0., 0., 13.5/27., 13.5/27., 0., 0., +// 0., 0., 0., 0., 6.75/27., 6.75/27., 6.75/27., 6.75/27., +// 0., 0., 0., 0., 13.5/27., 0., 0., 13.5/27., + }; + + static const double q1_into_q1_refined_5[] = + { +// 6.75/27., 6.75/27., 0., 0., 6.75/27., 6.75/27., 0., 0., +// 0., 13.5/27., 0., 0., 0., 13.5/27., 0., 0., +// 0., 6.75/27., 6.75/27., 0., 0., 6.75/27., 6.75/27., 0., +// 3.375/27., 3.375/27., 3.375/27., 3.375/27., 3.375/27., 3.375/27., 3.375/27., 3.375/27., +// 0., 0., 0., 0., 13.5/27., 13.5/27., 0., 0., +// 0., 0., 0., 0., 0., 1., 0., 0., +// 0., 0., 0., 0., 0., 13.5/27., 13.5/27., 0., +// 0., 0., 0., 0., 6.75/27., 6.75/27., 6.75/27., 6.75/27., + }; + + static const double q1_into_q1_refined_6[] = + { +// 3.375/27., 3.375/27., 3.375/27., 3.375/27., 3.375/27., 3.375/27., 3.375/27., 3.375/27., +// 0., 6.75/27., 6.75/27., 0., 0., 6.75/27., 6.75/27., 0., +// 0., 0., 13.5/27., 0., 0., 0., 13.5/27., 0., +// 0., 0., 6.75/27., 6.75/27., 0., 0., 6.75/27., 6.75/27., +// 0., 0., 0., 0., 6.75/27., 6.75/27., 6.75/27., 6.75/27., +// 0., 0., 0., 0., 0., 13.5/27., 13.5/27., 0., +// 0., 0., 0., 0., 0., 0., 1., 0., +// 0., 0., 0., 0., 0., 0., 13.5/27., 13.5/27., + }; + + static const double q1_into_q1_refined_7[] = + { +// 6.75/27., 0., 0., 6.75/27., 6.75/27., 0., 0., 6.75/27., +// 3.375/27., 3.375/27., 3.375/27., 3.375/27., 3.375/27., 3.375/27., 3.375/27., 3.375/27., +// 0., 0., 6.75/27., 6.75/27., 0., 0., 6.75/27., 6.75/27., +// 0., 0., 0., 13.5/27., 0., 0., 0., 13.5/27., +// 0., 0., 0., 0., 13.5/27., 0., 0., 13.5/27., +// 0., 0., 0., 0., 6.75/27., 6.75/27., 6.75/27., 6.75/27., +// 0., 0., 0., 0., 0., 0., 13.5/27., 13.5/27., +// 0., 0., 0., 0., 0., 0., 0., 1., + }; + +}; // namespace FE_Nedelec_3d + + +// embedding matrices + +template <> +const double * const +FE_Nedelec<3>::Matrices::embedding[][GeometryInfo<3>::children_per_cell] = +{ + { FE_Nedelec_3d::q1_into_q1_refined_0, FE_Nedelec_3d::q1_into_q1_refined_1, + FE_Nedelec_3d::q1_into_q1_refined_2, FE_Nedelec_3d::q1_into_q1_refined_3, + FE_Nedelec_3d::q1_into_q1_refined_4, FE_Nedelec_3d::q1_into_q1_refined_5, + FE_Nedelec_3d::q1_into_q1_refined_6, FE_Nedelec_3d::q1_into_q1_refined_7 } +}; + + +template <> +const unsigned int +FE_Nedelec<3>::Matrices::n_embedding_matrices + = sizeof(FE_Nedelec<3>::Matrices::embedding) / + sizeof(FE_Nedelec<3>::Matrices::embedding[0]); + + + +// Constraint matrices + +namespace FE_Nedelec_3d +{ + static const double constraint_q1[] = + { +//TODO[WB] +// .25,.25,.25,.25, +// .5,.5,0.,0., +// 0.,.5,.5,0., +// 0.,0.,.5,.5, +// .5,0.,0.,.5 + }; + +}; + + + +template <> +const double * const +FE_Nedelec<3>::Matrices::constraint_matrices[] = +{ + FE_Nedelec_3d::constraint_q1 +}; + + + +template <> +const unsigned int +FE_Nedelec<3>::Matrices::n_constraint_matrices + = sizeof(FE_Nedelec<3>::Matrices::constraint_matrices) / + sizeof(FE_Nedelec<3>::Matrices::constraint_matrices[0]); + + + +#else // #if deal_II_dimension +// On gcc2.95 on Alpha OSF1, the native assembler does not like empty +// files, so provide some dummy code +namespace { void dummy () {}; }; +#endif // #if deal_II_dimension == 3 diff --git a/deal.II/deal.II/source/fe/fe_q.cc b/deal.II/deal.II/source/fe/fe_q.cc index a483f7ba11..4f50b45931 100644 --- a/deal.II/deal.II/source/fe/fe_q.cc +++ b/deal.II/deal.II/source/fe/fe_q.cc @@ -23,29 +23,19 @@ #include - - template FE_Q::FE_Q (const unsigned int degree) : FiniteElement (FiniteElementData(get_dpo_vector(degree),1), - std::vector (1,false)), + std::vector (1,false), + std::vector >(FiniteElementData(get_dpo_vector(degree),1).dofs_per_cell, + std::vector(1,true))), degree(degree), renumber(dofs_per_cell, 0), renumber_inverse(dofs_per_cell, 0), face_renumber(dofs_per_face, 0), - poly(0) + polynomial_space(LagrangeEquidistant::generate_complete_basis(degree)) { - // Q0 elements cannot be - // continuous, use FE_DGQ(0) - // instead - Assert (degree>0, ExcNotImplemented()); - std::vector v; - for (unsigned int i=0;i<=degree;++i) - v.push_back(LagrangeEquidistant(degree,i)); - - poly = new TensorProductPolynomials (v); - // do some internal book-keeping on // cells and faces. if in 1d, the // face function is empty @@ -77,497 +67,553 @@ FE_Q::FE_Q (const unsigned int degree) // then fill restriction // matrices. they are hardcoded for - // the first few elements + // the first few elements. in + // contrast to the other matrices, + // these are not stored in the + // files fe_q_[123]d.cc, since they + // contain only a rather small + // number of zeros, and storing + // them element-wise is more + // expensive than just setting the + // nonzero elements as done here switch (dim) { - case 1: - switch (degree) - { - case 1: - restriction[0](0,0) = 1; - restriction[1](1,1) = 1; - break; - case 2: - restriction[0](0,0) = 1; - restriction[0](2,1) = 1; - restriction[1](1,1) = 1; - restriction[1](1,1) = 1; - break; - case 3: - restriction[0](0,0) = 1; - restriction[0](2,3) = 1; - restriction[1](1,1) = 1; - restriction[1](3,2) = 1; - break; - case 4: - restriction[0](0,0) = 1; - restriction[0](2,3) = 1; - restriction[0](3,1) = 1; - restriction[1](1,1) = 1; - restriction[1](3,0) = 1; - restriction[1](4,3) = 1; - break; - default: - for (unsigned int i=0; i::children_per_cell;++i) - restriction[i].reinit(0,0); - } - break; - case 2: - switch (degree) - { - case 1: - restriction[0](0,0) = 1; - restriction[1](1,1) = 1; - restriction[2](2,2) = 1; - restriction[3](3,3) = 1; - break; - case 2: - restriction[0](0,0) = 1; - restriction[0](4,1) = 1; - restriction[0](7,3) = 1; - restriction[0](8,2) = 1; - restriction[1](1,1) = 1; - restriction[1](4,0) = 1; - restriction[1](5,2) = 1; - restriction[1](8,3) = 1; - restriction[2](2,2) = 1; - restriction[2](5,1) = 1; - restriction[2](6,3) = 1; - restriction[2](8,0) = 1; - restriction[3](3,3) = 1; - restriction[3](6,2) = 1; - restriction[3](7,0) = 1; - restriction[3](8,1) = 1; - break; - case 3: - restriction[0](0,0) = 1; - restriction[0](4,5) = 1; - restriction[0](10,11) = 1; - restriction[0](12,15) = 1; - restriction[1](1,1) = 1; - restriction[1](5,4) = 1; - restriction[1](6,7) = 1; - restriction[1](13,14) = 1; - restriction[2](2,2) = 1; - restriction[2](7,6) = 1; - restriction[2](9,8) = 1; - restriction[2](15,12) = 1; - restriction[3](3,3) = 1; - restriction[3](8,9) = 1; - restriction[3](11,10) = 1; - restriction[3](14,13) = 1; - break; - case 4: - restriction[0](0,0) = 1; - restriction[0](4,5) = 1; - restriction[0](5,1) = 1; - restriction[0](13,14) = 1; - restriction[0](14,3) = 1; - restriction[0](16,20) = 1; - restriction[0](17,8) = 1; - restriction[0](19,11) = 1; - restriction[0](20,2) = 1; - restriction[1](1,1) = 1; - restriction[1](5,0) = 1; - restriction[1](6,5) = 1; - restriction[1](7,8) = 1; - restriction[1](8,2) = 1; - restriction[1](17,14) = 1; - restriction[1](18,20) = 1; - restriction[1](20,3) = 1; - restriction[1](21,11) = 1; - restriction[2](2,2) = 1; - restriction[2](8,1) = 1; - restriction[2](9,8) = 1; - restriction[2](11,3) = 1; - restriction[2](12,11) = 1; - restriction[2](20,0) = 1; - restriction[2](21,5) = 1; - restriction[2](23,14) = 1; - restriction[2](24,20) = 1; - restriction[3](3,3) = 1; - restriction[3](10,11) = 1; - restriction[3](11,2) = 1; - restriction[3](14,0) = 1; - restriction[3](15,14) = 1; - restriction[3](19,5) = 1; - restriction[3](20,1) = 1; - restriction[3](22,20) = 1; - restriction[3](23,8) = 1; - break; - default: - for (unsigned int i=0; i::children_per_cell;++i) - restriction[i].reinit(0,0); - } - break; - case 3: - switch (degree) - { - case 1: - restriction[0](0,0) = 1; - restriction[1](1,1) = 1; - restriction[2](2,2) = 1; - restriction[3](3,3) = 1; - restriction[4](4,4) = 1; - restriction[5](5,5) = 1; - restriction[6](6,6) = 1; - restriction[7](7,7) = 1; - break; - case 2: - restriction[0](0,0) = 1; - restriction[0](8,1) = 1; - restriction[0](11,3) = 1; - restriction[0](16,4) = 1; - restriction[0](20,2) = 1; - restriction[0](22,5) = 1; - restriction[0](25,7) = 1; - restriction[0](26,6) = 1; - restriction[1](1,1) = 1; - restriction[1](8,0) = 1; - restriction[1](9,2) = 1; - restriction[1](17,5) = 1; - restriction[1](20,3) = 1; - restriction[1](22,4) = 1; - restriction[1](23,6) = 1; - restriction[1](26,7) = 1; - restriction[2](2,2) = 1; - restriction[2](9,1) = 1; - restriction[2](10,3) = 1; - restriction[2](18,6) = 1; - restriction[2](20,0) = 1; - restriction[2](23,5) = 1; - restriction[2](24,7) = 1; - restriction[2](26,4) = 1; - restriction[3](3,3) = 1; - restriction[3](10,2) = 1; - restriction[3](11,0) = 1; - restriction[3](19,7) = 1; - restriction[3](20,1) = 1; - restriction[3](24,6) = 1; - restriction[3](25,4) = 1; - restriction[3](26,5) = 1; - restriction[4](4,4) = 1; - restriction[4](12,5) = 1; - restriction[4](15,7) = 1; - restriction[4](16,0) = 1; - restriction[4](21,6) = 1; - restriction[4](22,1) = 1; - restriction[4](25,3) = 1; - restriction[4](26,2) = 1; - restriction[5](5,5) = 1; - restriction[5](12,4) = 1; - restriction[5](13,6) = 1; - restriction[5](17,1) = 1; - restriction[5](21,7) = 1; - restriction[5](22,0) = 1; - restriction[5](23,2) = 1; - restriction[5](26,3) = 1; - restriction[6](6,6) = 1; - restriction[6](13,5) = 1; - restriction[6](14,7) = 1; - restriction[6](18,2) = 1; - restriction[6](21,4) = 1; - restriction[6](23,1) = 1; - restriction[6](24,3) = 1; - restriction[6](26,0) = 1; - restriction[7](7,7) = 1; - restriction[7](14,6) = 1; - restriction[7](15,4) = 1; - restriction[7](19,3) = 1; - restriction[7](21,5) = 1; - restriction[7](24,2) = 1; - restriction[7](25,0) = 1; - restriction[7](26,1) = 1; - break; - case 3: - restriction[0](0,0) = 1; - restriction[0](8,9) = 1; - restriction[0](14,15) = 1; - restriction[0](24,25) = 1; - restriction[0](32,35) = 1; - restriction[0](40,43) = 1; - restriction[0](52,55) = 1; - restriction[0](56,63) = 1; - restriction[1](1,1) = 1; - restriction[1](9,8) = 1; - restriction[1](10,11) = 1; - restriction[1](26,27) = 1; - restriction[1](33,34) = 1; - restriction[1](41,42) = 1; - restriction[1](44,47) = 1; - restriction[1](57,62) = 1; - restriction[2](2,2) = 1; - restriction[2](11,10) = 1; - restriction[2](13,12) = 1; - restriction[2](28,29) = 1; - restriction[2](35,32) = 1; - restriction[2](46,45) = 1; - restriction[2](49,50) = 1; - restriction[2](61,58) = 1; - restriction[3](3,3) = 1; - restriction[3](12,13) = 1; - restriction[3](15,14) = 1; - restriction[3](30,31) = 1; - restriction[3](34,33) = 1; - restriction[3](48,51) = 1; - restriction[3](54,53) = 1; - restriction[3](60,59) = 1; - restriction[4](4,4) = 1; - restriction[4](16,17) = 1; - restriction[4](22,23) = 1; - restriction[4](25,24) = 1; - restriction[4](36,39) = 1; - restriction[4](42,41) = 1; - restriction[4](53,54) = 1; - restriction[4](58,61) = 1; - restriction[5](5,5) = 1; - restriction[5](17,16) = 1; - restriction[5](18,19) = 1; - restriction[5](27,26) = 1; - restriction[5](37,38) = 1; - restriction[5](43,40) = 1; - restriction[5](45,46) = 1; - restriction[5](59,60) = 1; - restriction[6](6,6) = 1; - restriction[6](19,18) = 1; - restriction[6](21,20) = 1; - restriction[6](29,28) = 1; - restriction[6](39,36) = 1; - restriction[6](47,44) = 1; - restriction[6](51,48) = 1; - restriction[6](63,56) = 1; - restriction[7](7,7) = 1; - restriction[7](20,21) = 1; - restriction[7](23,22) = 1; - restriction[7](31,30) = 1; - restriction[7](38,37) = 1; - restriction[7](50,49) = 1; - restriction[7](55,52) = 1; - restriction[7](62,57) = 1; - break; - case 4: - restriction[0](0,0) = 1; - restriction[0](8,9) = 1; - restriction[0](9,1) = 1; - restriction[0](17,18) = 1; - restriction[0](18,3) = 1; - restriction[0](32,33) = 1; - restriction[0](33,4) = 1; - restriction[0](44,48) = 1; - restriction[0](45,12) = 1; - restriction[0](47,15) = 1; - restriction[0](48,2) = 1; - restriction[0](62,66) = 1; - restriction[0](63,36) = 1; - restriction[0](65,21) = 1; - restriction[0](66,5) = 1; - restriction[0](89,93) = 1; - restriction[0](90,30) = 1; - restriction[0](92,42) = 1; - restriction[0](93,7) = 1; - restriction[0](98,111) = 1; - restriction[0](99,75) = 1; - restriction[0](101,57) = 1; - restriction[0](102,24) = 1; - restriction[0](107,84) = 1; - restriction[0](108,39) = 1; - restriction[0](110,27) = 1; - restriction[0](111,6) = 1; - restriction[1](1,1) = 1; - restriction[1](9,0) = 1; - restriction[1](10,9) = 1; - restriction[1](11,12) = 1; - restriction[1](12,2) = 1; - restriction[1](35,36) = 1; - restriction[1](36,5) = 1; - restriction[1](45,18) = 1; - restriction[1](46,48) = 1; - restriction[1](48,3) = 1; - restriction[1](49,15) = 1; - restriction[1](63,33) = 1; - restriction[1](64,66) = 1; - restriction[1](66,4) = 1; - restriction[1](67,21) = 1; - restriction[1](71,75) = 1; - restriction[1](72,24) = 1; - restriction[1](74,39) = 1; - restriction[1](75,6) = 1; - restriction[1](99,93) = 1; - restriction[1](100,111) = 1; - restriction[1](102,30) = 1; - restriction[1](103,57) = 1; - restriction[1](108,42) = 1; - restriction[1](109,84) = 1; - restriction[1](111,7) = 1; - restriction[1](112,27) = 1; - restriction[2](2,2) = 1; - restriction[2](12,1) = 1; - restriction[2](13,12) = 1; - restriction[2](15,3) = 1; - restriction[2](16,15) = 1; - restriction[2](38,39) = 1; - restriction[2](39,6) = 1; - restriction[2](48,0) = 1; - restriction[2](49,9) = 1; - restriction[2](51,18) = 1; - restriction[2](52,48) = 1; - restriction[2](74,36) = 1; - restriction[2](75,5) = 1; - restriction[2](77,75) = 1; - restriction[2](78,24) = 1; - restriction[2](81,42) = 1; - restriction[2](82,84) = 1; - restriction[2](84,7) = 1; - restriction[2](85,27) = 1; - restriction[2](108,33) = 1; - restriction[2](109,66) = 1; - restriction[2](111,4) = 1; - restriction[2](112,21) = 1; - restriction[2](117,93) = 1; - restriction[2](118,111) = 1; - restriction[2](120,30) = 1; - restriction[2](121,57) = 1; - restriction[3](3,3) = 1; - restriction[3](14,15) = 1; - restriction[3](15,2) = 1; - restriction[3](18,0) = 1; - restriction[3](19,18) = 1; - restriction[3](41,42) = 1; - restriction[3](42,7) = 1; - restriction[3](47,9) = 1; - restriction[3](48,1) = 1; - restriction[3](50,48) = 1; - restriction[3](51,12) = 1; - restriction[3](80,84) = 1; - restriction[3](81,39) = 1; - restriction[3](83,27) = 1; - restriction[3](84,6) = 1; - restriction[3](92,33) = 1; - restriction[3](93,4) = 1; - restriction[3](95,93) = 1; - restriction[3](96,30) = 1; - restriction[3](107,66) = 1; - restriction[3](108,36) = 1; - restriction[3](110,21) = 1; - restriction[3](111,5) = 1; - restriction[3](116,111) = 1; - restriction[3](117,75) = 1; - restriction[3](119,57) = 1; - restriction[3](120,24) = 1; - restriction[4](4,4) = 1; - restriction[4](20,21) = 1; - restriction[4](21,5) = 1; - restriction[4](29,30) = 1; - restriction[4](30,7) = 1; - restriction[4](33,0) = 1; - restriction[4](34,33) = 1; - restriction[4](53,57) = 1; - restriction[4](54,24) = 1; - restriction[4](56,27) = 1; - restriction[4](57,6) = 1; - restriction[4](65,9) = 1; - restriction[4](66,1) = 1; - restriction[4](68,66) = 1; - restriction[4](69,36) = 1; - restriction[4](90,18) = 1; - restriction[4](91,93) = 1; - restriction[4](93,3) = 1; - restriction[4](94,42) = 1; - restriction[4](101,48) = 1; - restriction[4](102,12) = 1; - restriction[4](104,111) = 1; - restriction[4](105,75) = 1; - restriction[4](110,15) = 1; - restriction[4](111,2) = 1; - restriction[4](113,84) = 1; - restriction[4](114,39) = 1; - restriction[5](5,5) = 1; - restriction[5](21,4) = 1; - restriction[5](22,21) = 1; - restriction[5](23,24) = 1; - restriction[5](24,6) = 1; - restriction[5](36,1) = 1; - restriction[5](37,36) = 1; - restriction[5](54,30) = 1; - restriction[5](55,57) = 1; - restriction[5](57,7) = 1; - restriction[5](58,27) = 1; - restriction[5](66,0) = 1; - restriction[5](67,9) = 1; - restriction[5](69,33) = 1; - restriction[5](70,66) = 1; - restriction[5](72,12) = 1; - restriction[5](73,75) = 1; - restriction[5](75,2) = 1; - restriction[5](76,39) = 1; - restriction[5](102,18) = 1; - restriction[5](103,48) = 1; - restriction[5](105,93) = 1; - restriction[5](106,111) = 1; - restriction[5](111,3) = 1; - restriction[5](112,15) = 1; - restriction[5](114,42) = 1; - restriction[5](115,84) = 1; - restriction[6](6,6) = 1; - restriction[6](24,5) = 1; - restriction[6](25,24) = 1; - restriction[6](27,7) = 1; - restriction[6](28,27) = 1; - restriction[6](39,2) = 1; - restriction[6](40,39) = 1; - restriction[6](57,4) = 1; - restriction[6](58,21) = 1; - restriction[6](60,30) = 1; - restriction[6](61,57) = 1; - restriction[6](75,1) = 1; - restriction[6](76,36) = 1; - restriction[6](78,12) = 1; - restriction[6](79,75) = 1; - restriction[6](84,3) = 1; - restriction[6](85,15) = 1; - restriction[6](87,42) = 1; - restriction[6](88,84) = 1; - restriction[6](111,0) = 1; - restriction[6](112,9) = 1; - restriction[6](114,33) = 1; - restriction[6](115,66) = 1; - restriction[6](120,18) = 1; - restriction[6](121,48) = 1; - restriction[6](123,93) = 1; - restriction[6](124,111) = 1; - restriction[7](7,7) = 1; - restriction[7](26,27) = 1; - restriction[7](27,6) = 1; - restriction[7](30,4) = 1; - restriction[7](31,30) = 1; - restriction[7](42,3) = 1; - restriction[7](43,42) = 1; - restriction[7](56,21) = 1; - restriction[7](57,5) = 1; - restriction[7](59,57) = 1; - restriction[7](60,24) = 1; - restriction[7](83,15) = 1; - restriction[7](84,2) = 1; - restriction[7](86,84) = 1; - restriction[7](87,39) = 1; - restriction[7](93,0) = 1; - restriction[7](94,33) = 1; - restriction[7](96,18) = 1; - restriction[7](97,93) = 1; - restriction[7](110,9) = 1; - restriction[7](111,1) = 1; - restriction[7](113,66) = 1; - restriction[7](114,36) = 1; - restriction[7](119,48) = 1; - restriction[7](120,12) = 1; - restriction[7](122,111) = 1; - restriction[7](123,75) = 1; - break; - default: - for (unsigned int i=0; i::children_per_cell;++i) - restriction[i].reinit(0,0); - } - break; + case 1: // 1d + { + switch (degree) + { + case 1: + restriction[0](0,0) = 1; + restriction[1](1,1) = 1; + break; + case 2: + restriction[0](0,0) = 1; + restriction[0](2,1) = 1; + restriction[1](1,1) = 1; + restriction[1](1,1) = 1; + break; + case 3: + restriction[0](0,0) = 1; + restriction[0](2,3) = 1; + restriction[1](1,1) = 1; + restriction[1](3,2) = 1; + break; + case 4: + restriction[0](0,0) = 1; + restriction[0](2,3) = 1; + restriction[0](3,1) = 1; + restriction[1](1,1) = 1; + restriction[1](3,0) = 1; + restriction[1](4,3) = 1; + break; + { + // in case we don't + // have the matrices + // (yet), set them to + // impossible + // values. this does + // not prevent the use + // of this FE, but will + // prevent the use of + // these matrices + for (unsigned int i=0; + i::children_per_cell; + ++i) + restriction[i].reinit(0,0); + }; + } + break; + }; + + case 2: // 2d + { + switch (degree) + { + case 1: + restriction[0](0,0) = 1; + restriction[1](1,1) = 1; + restriction[2](2,2) = 1; + restriction[3](3,3) = 1; + break; + case 2: + restriction[0](0,0) = 1; + restriction[0](4,1) = 1; + restriction[0](7,3) = 1; + restriction[0](8,2) = 1; + restriction[1](1,1) = 1; + restriction[1](4,0) = 1; + restriction[1](5,2) = 1; + restriction[1](8,3) = 1; + restriction[2](2,2) = 1; + restriction[2](5,1) = 1; + restriction[2](6,3) = 1; + restriction[2](8,0) = 1; + restriction[3](3,3) = 1; + restriction[3](6,2) = 1; + restriction[3](7,0) = 1; + restriction[3](8,1) = 1; + break; + case 3: + restriction[0](0,0) = 1; + restriction[0](4,5) = 1; + restriction[0](10,11) = 1; + restriction[0](12,15) = 1; + restriction[1](1,1) = 1; + restriction[1](5,4) = 1; + restriction[1](6,7) = 1; + restriction[1](13,14) = 1; + restriction[2](2,2) = 1; + restriction[2](7,6) = 1; + restriction[2](9,8) = 1; + restriction[2](15,12) = 1; + restriction[3](3,3) = 1; + restriction[3](8,9) = 1; + restriction[3](11,10) = 1; + restriction[3](14,13) = 1; + break; + case 4: + restriction[0](0,0) = 1; + restriction[0](4,5) = 1; + restriction[0](5,1) = 1; + restriction[0](13,14) = 1; + restriction[0](14,3) = 1; + restriction[0](16,20) = 1; + restriction[0](17,8) = 1; + restriction[0](19,11) = 1; + restriction[0](20,2) = 1; + restriction[1](1,1) = 1; + restriction[1](5,0) = 1; + restriction[1](6,5) = 1; + restriction[1](7,8) = 1; + restriction[1](8,2) = 1; + restriction[1](17,14) = 1; + restriction[1](18,20) = 1; + restriction[1](20,3) = 1; + restriction[1](21,11) = 1; + restriction[2](2,2) = 1; + restriction[2](8,1) = 1; + restriction[2](9,8) = 1; + restriction[2](11,3) = 1; + restriction[2](12,11) = 1; + restriction[2](20,0) = 1; + restriction[2](21,5) = 1; + restriction[2](23,14) = 1; + restriction[2](24,20) = 1; + restriction[3](3,3) = 1; + restriction[3](10,11) = 1; + restriction[3](11,2) = 1; + restriction[3](14,0) = 1; + restriction[3](15,14) = 1; + restriction[3](19,5) = 1; + restriction[3](20,1) = 1; + restriction[3](22,20) = 1; + restriction[3](23,8) = 1; + break; + + default: + { + // in case we don't + // have the matrices + // (yet), set them to + // impossible + // values. this does + // not prevent the use + // of this FE, but will + // prevent the use of + // these matrices + for (unsigned int i=0; + i::children_per_cell; + ++i) + restriction[i].reinit(0,0); + }; + } + break; + }; + + case 3: // 3d + { + switch (degree) + { + case 1: + restriction[0](0,0) = 1; + restriction[1](1,1) = 1; + restriction[2](2,2) = 1; + restriction[3](3,3) = 1; + restriction[4](4,4) = 1; + restriction[5](5,5) = 1; + restriction[6](6,6) = 1; + restriction[7](7,7) = 1; + break; + case 2: + restriction[0](0,0) = 1; + restriction[0](8,1) = 1; + restriction[0](11,3) = 1; + restriction[0](16,4) = 1; + restriction[0](20,2) = 1; + restriction[0](22,5) = 1; + restriction[0](25,7) = 1; + restriction[0](26,6) = 1; + restriction[1](1,1) = 1; + restriction[1](8,0) = 1; + restriction[1](9,2) = 1; + restriction[1](17,5) = 1; + restriction[1](20,3) = 1; + restriction[1](22,4) = 1; + restriction[1](23,6) = 1; + restriction[1](26,7) = 1; + restriction[2](2,2) = 1; + restriction[2](9,1) = 1; + restriction[2](10,3) = 1; + restriction[2](18,6) = 1; + restriction[2](20,0) = 1; + restriction[2](23,5) = 1; + restriction[2](24,7) = 1; + restriction[2](26,4) = 1; + restriction[3](3,3) = 1; + restriction[3](10,2) = 1; + restriction[3](11,0) = 1; + restriction[3](19,7) = 1; + restriction[3](20,1) = 1; + restriction[3](24,6) = 1; + restriction[3](25,4) = 1; + restriction[3](26,5) = 1; + restriction[4](4,4) = 1; + restriction[4](12,5) = 1; + restriction[4](15,7) = 1; + restriction[4](16,0) = 1; + restriction[4](21,6) = 1; + restriction[4](22,1) = 1; + restriction[4](25,3) = 1; + restriction[4](26,2) = 1; + restriction[5](5,5) = 1; + restriction[5](12,4) = 1; + restriction[5](13,6) = 1; + restriction[5](17,1) = 1; + restriction[5](21,7) = 1; + restriction[5](22,0) = 1; + restriction[5](23,2) = 1; + restriction[5](26,3) = 1; + restriction[6](6,6) = 1; + restriction[6](13,5) = 1; + restriction[6](14,7) = 1; + restriction[6](18,2) = 1; + restriction[6](21,4) = 1; + restriction[6](23,1) = 1; + restriction[6](24,3) = 1; + restriction[6](26,0) = 1; + restriction[7](7,7) = 1; + restriction[7](14,6) = 1; + restriction[7](15,4) = 1; + restriction[7](19,3) = 1; + restriction[7](21,5) = 1; + restriction[7](24,2) = 1; + restriction[7](25,0) = 1; + restriction[7](26,1) = 1; + break; + case 3: + restriction[0](0,0) = 1; + restriction[0](8,9) = 1; + restriction[0](14,15) = 1; + restriction[0](24,25) = 1; + restriction[0](32,35) = 1; + restriction[0](40,43) = 1; + restriction[0](52,55) = 1; + restriction[0](56,63) = 1; + restriction[1](1,1) = 1; + restriction[1](9,8) = 1; + restriction[1](10,11) = 1; + restriction[1](26,27) = 1; + restriction[1](33,34) = 1; + restriction[1](41,42) = 1; + restriction[1](44,47) = 1; + restriction[1](57,62) = 1; + restriction[2](2,2) = 1; + restriction[2](11,10) = 1; + restriction[2](13,12) = 1; + restriction[2](28,29) = 1; + restriction[2](35,32) = 1; + restriction[2](46,45) = 1; + restriction[2](49,50) = 1; + restriction[2](61,58) = 1; + restriction[3](3,3) = 1; + restriction[3](12,13) = 1; + restriction[3](15,14) = 1; + restriction[3](30,31) = 1; + restriction[3](34,33) = 1; + restriction[3](48,51) = 1; + restriction[3](54,53) = 1; + restriction[3](60,59) = 1; + restriction[4](4,4) = 1; + restriction[4](16,17) = 1; + restriction[4](22,23) = 1; + restriction[4](25,24) = 1; + restriction[4](36,39) = 1; + restriction[4](42,41) = 1; + restriction[4](53,54) = 1; + restriction[4](58,61) = 1; + restriction[5](5,5) = 1; + restriction[5](17,16) = 1; + restriction[5](18,19) = 1; + restriction[5](27,26) = 1; + restriction[5](37,38) = 1; + restriction[5](43,40) = 1; + restriction[5](45,46) = 1; + restriction[5](59,60) = 1; + restriction[6](6,6) = 1; + restriction[6](19,18) = 1; + restriction[6](21,20) = 1; + restriction[6](29,28) = 1; + restriction[6](39,36) = 1; + restriction[6](47,44) = 1; + restriction[6](51,48) = 1; + restriction[6](63,56) = 1; + restriction[7](7,7) = 1; + restriction[7](20,21) = 1; + restriction[7](23,22) = 1; + restriction[7](31,30) = 1; + restriction[7](38,37) = 1; + restriction[7](50,49) = 1; + restriction[7](55,52) = 1; + restriction[7](62,57) = 1; + break; + case 4: + restriction[0](0,0) = 1; + restriction[0](8,9) = 1; + restriction[0](9,1) = 1; + restriction[0](17,18) = 1; + restriction[0](18,3) = 1; + restriction[0](32,33) = 1; + restriction[0](33,4) = 1; + restriction[0](44,48) = 1; + restriction[0](45,12) = 1; + restriction[0](47,15) = 1; + restriction[0](48,2) = 1; + restriction[0](62,66) = 1; + restriction[0](63,36) = 1; + restriction[0](65,21) = 1; + restriction[0](66,5) = 1; + restriction[0](89,93) = 1; + restriction[0](90,30) = 1; + restriction[0](92,42) = 1; + restriction[0](93,7) = 1; + restriction[0](98,111) = 1; + restriction[0](99,75) = 1; + restriction[0](101,57) = 1; + restriction[0](102,24) = 1; + restriction[0](107,84) = 1; + restriction[0](108,39) = 1; + restriction[0](110,27) = 1; + restriction[0](111,6) = 1; + restriction[1](1,1) = 1; + restriction[1](9,0) = 1; + restriction[1](10,9) = 1; + restriction[1](11,12) = 1; + restriction[1](12,2) = 1; + restriction[1](35,36) = 1; + restriction[1](36,5) = 1; + restriction[1](45,18) = 1; + restriction[1](46,48) = 1; + restriction[1](48,3) = 1; + restriction[1](49,15) = 1; + restriction[1](63,33) = 1; + restriction[1](64,66) = 1; + restriction[1](66,4) = 1; + restriction[1](67,21) = 1; + restriction[1](71,75) = 1; + restriction[1](72,24) = 1; + restriction[1](74,39) = 1; + restriction[1](75,6) = 1; + restriction[1](99,93) = 1; + restriction[1](100,111) = 1; + restriction[1](102,30) = 1; + restriction[1](103,57) = 1; + restriction[1](108,42) = 1; + restriction[1](109,84) = 1; + restriction[1](111,7) = 1; + restriction[1](112,27) = 1; + restriction[2](2,2) = 1; + restriction[2](12,1) = 1; + restriction[2](13,12) = 1; + restriction[2](15,3) = 1; + restriction[2](16,15) = 1; + restriction[2](38,39) = 1; + restriction[2](39,6) = 1; + restriction[2](48,0) = 1; + restriction[2](49,9) = 1; + restriction[2](51,18) = 1; + restriction[2](52,48) = 1; + restriction[2](74,36) = 1; + restriction[2](75,5) = 1; + restriction[2](77,75) = 1; + restriction[2](78,24) = 1; + restriction[2](81,42) = 1; + restriction[2](82,84) = 1; + restriction[2](84,7) = 1; + restriction[2](85,27) = 1; + restriction[2](108,33) = 1; + restriction[2](109,66) = 1; + restriction[2](111,4) = 1; + restriction[2](112,21) = 1; + restriction[2](117,93) = 1; + restriction[2](118,111) = 1; + restriction[2](120,30) = 1; + restriction[2](121,57) = 1; + restriction[3](3,3) = 1; + restriction[3](14,15) = 1; + restriction[3](15,2) = 1; + restriction[3](18,0) = 1; + restriction[3](19,18) = 1; + restriction[3](41,42) = 1; + restriction[3](42,7) = 1; + restriction[3](47,9) = 1; + restriction[3](48,1) = 1; + restriction[3](50,48) = 1; + restriction[3](51,12) = 1; + restriction[3](80,84) = 1; + restriction[3](81,39) = 1; + restriction[3](83,27) = 1; + restriction[3](84,6) = 1; + restriction[3](92,33) = 1; + restriction[3](93,4) = 1; + restriction[3](95,93) = 1; + restriction[3](96,30) = 1; + restriction[3](107,66) = 1; + restriction[3](108,36) = 1; + restriction[3](110,21) = 1; + restriction[3](111,5) = 1; + restriction[3](116,111) = 1; + restriction[3](117,75) = 1; + restriction[3](119,57) = 1; + restriction[3](120,24) = 1; + restriction[4](4,4) = 1; + restriction[4](20,21) = 1; + restriction[4](21,5) = 1; + restriction[4](29,30) = 1; + restriction[4](30,7) = 1; + restriction[4](33,0) = 1; + restriction[4](34,33) = 1; + restriction[4](53,57) = 1; + restriction[4](54,24) = 1; + restriction[4](56,27) = 1; + restriction[4](57,6) = 1; + restriction[4](65,9) = 1; + restriction[4](66,1) = 1; + restriction[4](68,66) = 1; + restriction[4](69,36) = 1; + restriction[4](90,18) = 1; + restriction[4](91,93) = 1; + restriction[4](93,3) = 1; + restriction[4](94,42) = 1; + restriction[4](101,48) = 1; + restriction[4](102,12) = 1; + restriction[4](104,111) = 1; + restriction[4](105,75) = 1; + restriction[4](110,15) = 1; + restriction[4](111,2) = 1; + restriction[4](113,84) = 1; + restriction[4](114,39) = 1; + restriction[5](5,5) = 1; + restriction[5](21,4) = 1; + restriction[5](22,21) = 1; + restriction[5](23,24) = 1; + restriction[5](24,6) = 1; + restriction[5](36,1) = 1; + restriction[5](37,36) = 1; + restriction[5](54,30) = 1; + restriction[5](55,57) = 1; + restriction[5](57,7) = 1; + restriction[5](58,27) = 1; + restriction[5](66,0) = 1; + restriction[5](67,9) = 1; + restriction[5](69,33) = 1; + restriction[5](70,66) = 1; + restriction[5](72,12) = 1; + restriction[5](73,75) = 1; + restriction[5](75,2) = 1; + restriction[5](76,39) = 1; + restriction[5](102,18) = 1; + restriction[5](103,48) = 1; + restriction[5](105,93) = 1; + restriction[5](106,111) = 1; + restriction[5](111,3) = 1; + restriction[5](112,15) = 1; + restriction[5](114,42) = 1; + restriction[5](115,84) = 1; + restriction[6](6,6) = 1; + restriction[6](24,5) = 1; + restriction[6](25,24) = 1; + restriction[6](27,7) = 1; + restriction[6](28,27) = 1; + restriction[6](39,2) = 1; + restriction[6](40,39) = 1; + restriction[6](57,4) = 1; + restriction[6](58,21) = 1; + restriction[6](60,30) = 1; + restriction[6](61,57) = 1; + restriction[6](75,1) = 1; + restriction[6](76,36) = 1; + restriction[6](78,12) = 1; + restriction[6](79,75) = 1; + restriction[6](84,3) = 1; + restriction[6](85,15) = 1; + restriction[6](87,42) = 1; + restriction[6](88,84) = 1; + restriction[6](111,0) = 1; + restriction[6](112,9) = 1; + restriction[6](114,33) = 1; + restriction[6](115,66) = 1; + restriction[6](120,18) = 1; + restriction[6](121,48) = 1; + restriction[6](123,93) = 1; + restriction[6](124,111) = 1; + restriction[7](7,7) = 1; + restriction[7](26,27) = 1; + restriction[7](27,6) = 1; + restriction[7](30,4) = 1; + restriction[7](31,30) = 1; + restriction[7](42,3) = 1; + restriction[7](43,42) = 1; + restriction[7](56,21) = 1; + restriction[7](57,5) = 1; + restriction[7](59,57) = 1; + restriction[7](60,24) = 1; + restriction[7](83,15) = 1; + restriction[7](84,2) = 1; + restriction[7](86,84) = 1; + restriction[7](87,39) = 1; + restriction[7](93,0) = 1; + restriction[7](94,33) = 1; + restriction[7](96,18) = 1; + restriction[7](97,93) = 1; + restriction[7](110,9) = 1; + restriction[7](111,1) = 1; + restriction[7](113,66) = 1; + restriction[7](114,36) = 1; + restriction[7](119,48) = 1; + restriction[7](120,12) = 1; + restriction[7](122,111) = 1; + restriction[7](123,75) = 1; + break; + default: + { + // in case we don't + // have the matrices + // (yet), set them to + // impossible + // values. this does + // not prevent the use + // of this FE, but will + // prevent the use of + // these matrices + for (unsigned int i=0; + i::children_per_cell; + ++i) + restriction[i].reinit(0,0); + }; + } + break; + }; + default: Assert (false,ExcNotImplemented()); } @@ -580,17 +626,6 @@ FE_Q::FE_Q (const unsigned int degree) -template -FE_Q::~FE_Q () -{ - // delete poly member and set it to - // zero to prevent accidental use - delete poly; - poly = 0; -} - - - template FiniteElement * FE_Q::clone() const @@ -605,7 +640,20 @@ double FE_Q::shape_value (const unsigned int i, const Point &p) const { - return poly->compute_value(renumber_inverse[i], p); + Assert (i +double +FE_Q::shape_value_component (const unsigned int i, + const Point &p, + const unsigned int component) const +{ + Assert (i FE_Q::shape_grad (const unsigned int i, const Point &p) const { - return poly->compute_grad(renumber_inverse[i], p); + Assert (i +Tensor<1,dim> +FE_Q::shape_grad_component (const unsigned int i, + const Point &p, + const unsigned int component) const +{ + Assert (i FE_Q::shape_grad_grad (const unsigned int i, const Point &p) const { - return poly->compute_grad_grad(renumber_inverse[i], p); + Assert (i +Tensor<2,dim> +FE_Q::shape_grad_grad_component (const unsigned int i, + const Point &p, + const unsigned int component) const +{ + Assert (i UpdateFlags FE_Q::update_once (const UpdateFlags flags) const { - UpdateFlags out = update_default; - - if (flags & update_values) - out |= update_values; - - return out; + // for this kind of elements, only + // the values can be precomputed + // once and for all. set this flag + // if the values are requested at + // all + return (update_default | (flags & update_values)); } @@ -1006,10 +1082,9 @@ FE_Q::get_data (const UpdateFlags update_flags, const Mapping &mapping, const Quadrature &quadrature) const { + // generate a new data object and + // initialize some fields InternalData* data = new InternalData; - std::vector values(0); - std::vector > grads(0); - std::vector > grad_grads(0); // check what needs to be // initialized only once and what @@ -1021,7 +1096,15 @@ FE_Q::get_data (const UpdateFlags update_flags, const UpdateFlags flags(data->update_flags); const unsigned int n_q_points = quadrature.n_quadrature_points; - + + // some scratch arrays + std::vector values(0); + std::vector > grads(0); + std::vector > grad_grads(0); + + // initialize fields only if really + // necessary. otherwise, don't + // allocate memory if (flags & update_values) { values.resize (dofs_per_cell); @@ -1042,18 +1125,27 @@ FE_Q::get_data (const UpdateFlags update_flags, // that if (flags & update_second_derivatives) data->initialize_2nd (this, mapping, quadrature); - + + // next already fill those fields + // of which we have information by + // now. note that the shape + // gradients are only those on the + // unit cell, and need to be + // transformed when visiting an + // actual cell if (flags & (update_values | update_gradients)) for (unsigned int i=0; icompute(quadrature.point(i), values, grads, grad_grads); - for (unsigned int k=0; kshape_values[renumber[k]][i] = values[k]; - if (flags & update_gradients) - data->shape_gradients[renumber[k]][i] = grads[k]; - } + polynomial_space.compute(quadrature.point(i), + values, grads, grad_grads); + + if (flags & update_values) + for (unsigned int k=0; kshape_values[renumber[k]][i] = values[k]; + + if (flags & update_gradients) + for (unsigned int k=0; kshape_gradients[renumber[k]][i] = grads[k]; } return data; } @@ -1084,8 +1176,8 @@ FE_Q::fill_fe_values (const Mapping &mapping, for (unsigned int k=0; k::has_support_on_face (const unsigned int shape_index_, const unsigned int face_index) const { + Assert (shape_index_ < dofs_per_cell, + ExcIndexRange (shape_index_, 0, dofs_per_cell)); + Assert (face_index < GeometryInfo::faces_per_cell, + ExcIndexRange (face_index, 0, GeometryInfo::faces_per_cell)); + unsigned int shape_index = shape_index_; if (dim==1) diff --git a/deal.II/deal.II/source/fe/fe_system.cc b/deal.II/deal.II/source/fe/fe_system.cc index 0416090160..76614cb894 100644 --- a/deal.II/deal.II/source/fe/fe_system.cc +++ b/deal.II/deal.II/source/fe/fe_system.cc @@ -135,9 +135,11 @@ const unsigned int FESystem::invalid_face_number; template -FESystem::FESystem (const FiniteElement &fe, const unsigned int n_elements) : +FESystem::FESystem (const FiniteElement &fe, + const unsigned int n_elements) : FiniteElement (multiply_dof_numbers(fe, n_elements), - compute_restriction_is_additive_flags (fe, n_elements)), + compute_restriction_is_additive_flags (fe, n_elements), + compute_nonzero_components(fe, n_elements)), base_elements(1) { base_elements[0] = ElementPair(fe.clone(), n_elements); @@ -148,11 +150,15 @@ FESystem::FESystem (const FiniteElement &fe, const unsigned int n_elem template -FESystem::FESystem (const FiniteElement &fe1, const unsigned int n1, - const FiniteElement &fe2, const unsigned int n2) : +FESystem::FESystem (const FiniteElement &fe1, + const unsigned int n1, + const FiniteElement &fe2, + const unsigned int n2) : FiniteElement (multiply_dof_numbers(fe1, n1, fe2, n2), compute_restriction_is_additive_flags (fe1, n1, - fe2, n2)), + fe2, n2), + compute_nonzero_components(fe1, n1, + fe2, n2)), base_elements(2) { base_elements[0] = ElementPair(fe1.clone(), n1); @@ -165,15 +171,21 @@ FESystem::FESystem (const FiniteElement &fe1, const unsigned int n1, template -FESystem::FESystem (const FiniteElement &fe1, const unsigned int n1, - const FiniteElement &fe2, const unsigned int n2, - const FiniteElement &fe3, const unsigned int n3) : +FESystem::FESystem (const FiniteElement &fe1, + const unsigned int n1, + const FiniteElement &fe2, + const unsigned int n2, + const FiniteElement &fe3, + const unsigned int n3) : FiniteElement (multiply_dof_numbers(fe1, n1, fe2, n2, fe3, n3), compute_restriction_is_additive_flags (fe1, n1, fe2, n2, - fe3, n3)), + fe3, n3), + compute_nonzero_components(fe1, n1, + fe2, n2, + fe3, n3)), base_elements(3) { base_elements[0] = ElementPair(fe1.clone(), n1); @@ -236,12 +248,44 @@ double FESystem::shape_value (const unsigned int i, const Point &p) const { - Assert((i::ExcShapeFunctionNotPrimitive(i)); - std::pair comp = system_to_component_index(i); + return (base_element(system_to_base_table[i].first.first) + .shape_value(system_to_base_table[i].second, p)); +} + + + +template +double +FESystem::shape_value_component (const unsigned int i, + const Point &p, + const unsigned int component) const +{ + Assert (i FESystem::shape_grad (const unsigned int i, const Point &p) const { - Assert((i::ExcShapeFunctionNotPrimitive(i)); - std::pair comp = system_to_component_index(i); + return (base_element(system_to_base_table[i].first.first) + .shape_grad(system_to_base_table[i].second, p)); +} + + + +template +Tensor<1,dim> +FESystem::shape_grad_component (const unsigned int i, + const Point &p, + const unsigned int component) const +{ + Assert (i FESystem::shape_grad_grad (const unsigned int i, const Point &p) const { - Assert((i::ExcShapeFunctionNotPrimitive(i)); + + return (base_element(system_to_base_table[i].first.first) + .shape_grad_grad(system_to_base_table[i].second, p)); +} + - std::pair comp = system_to_component_index(i); +template +Tensor<2,dim> +FESystem::shape_grad_grad_component (const unsigned int i, + const Point &p, + const unsigned int component) const +{ + Assert (i::get_data (const UpdateFlags flags_, UpdateFlags sub_flags = flags; // if second derivatives through - // finite differencing is required, + // finite differencing are required, // then initialize some objects for // that data->compute_second_derivatives = flags & update_second_derivatives; @@ -398,7 +505,8 @@ FESystem::get_data (const UpdateFlags flags_, FEValuesData *base_data = new FEValuesData(); data->set_fe_values_data(base_no, base_data); } - data->update_flags=data->update_once | data->update_each; + data->update_flags = data->update_once | + data->update_each; return data; } @@ -509,7 +617,7 @@ FESystem::compute_fill (const Mapping &mapping, FEValuesData &base_data=fe_data.get_fe_values_data(base_no); const FiniteElement &base_fe=base_element(base_no); base_data.initialize(quadrature.n_quadrature_points, - base_fe.dofs_per_cell, + base_fe, base_update_flags); } } @@ -657,276 +765,462 @@ FESystem::compute_fill (const Mapping &mapping, template void -FESystem::build_cell_table() +FESystem::build_cell_tables() { unsigned total_index = 0; - for (unsigned int base=0 ; base < n_base_elements() ; ++base) + for (unsigned int base=0; base < n_base_elements(); ++base) for (unsigned int m = 0; m < element_multiplicity(base); ++m) for (unsigned int k=0; k + non_primitive_index (static_cast(-1), + static_cast(-1)); - // 1. Vertices + // First enumerate vertex indices, + // where we first enumerate all + // indices on the first vertex in + // the order of the base elements, + // then of the second vertex, etc total_index = 0; - for (unsigned int vertex_number= 0 ; - vertex_number < GeometryInfo::vertices_per_cell ; + for (unsigned int vertex_number=0; + vertex_number::vertices_per_cell; ++vertex_number) { unsigned comp_start = 0; - for(unsigned int base = 0; base < n_base_elements() ; - ++base) - { - for (unsigned int m = 0; m < element_multiplicity(base); ++m) + for(unsigned int base=0; base::lines_per_cell) && - (GeometryInfo::lines_per_cell > 0)); - ++line_number) - { - unsigned comp_start = 0; - for(unsigned int base = 0; base < n_base_elements() ; - ++base) - { - for (unsigned int m = 0; m < element_multiplicity(base); ++m) - { - for (unsigned int local_index = 0 ; - local_index < base_element(base).dofs_per_line ; - ++local_index) - { - system_to_component_table[total_index++] - = std::pair - (comp_start+m, - line_number*base_element(base).dofs_per_line - +local_index+base_element(base).first_line_index); - } - } - comp_start += element_multiplicity(base); - } - } + if (GeometryInfo::lines_per_cell > 0) + for (unsigned int line_number= 0; + line_number != GeometryInfo::lines_per_cell; + ++line_number) + { + unsigned comp_start = 0; + for (unsigned int base=0; base::quads_per_cell) && - (GeometryInfo::quads_per_cell > 0)); - ++quad_number) - { - unsigned int comp_start = 0; - for(unsigned int base = 0; base < n_base_elements() ; - ++base) - { - for (unsigned int m = 0; m < element_multiplicity(base); ++m) - { - for (unsigned int local_index = 0 ; - local_index < base_element(base).dofs_per_quad ; - ++local_index) - { - system_to_component_table[total_index++] - = std::make_pair (comp_start+m, - quad_number*base_element(base).dofs_per_quad - +local_index+base_element(base).first_quad_index); - } - } - comp_start += element_multiplicity(base); - } - } + if (GeometryInfo::quads_per_cell > 0) + for (unsigned int quad_number= 0; + quad_number != GeometryInfo::quads_per_cell; + ++quad_number) + { + unsigned int comp_start = 0; + for (unsigned int base=0; base::hexes_per_cell) && - (GeometryInfo::hexes_per_cell > 0)); - ++hex_number) + // 4. Hexes + if (GeometryInfo::hexes_per_cell > 0) + for (unsigned int hex_number= 0; + hex_number != GeometryInfo::hexes_per_cell; + ++hex_number) + { + unsigned int comp_start = 0; + for(unsigned int base=0; base dofs_per_component (n_components(), 0); + for (unsigned int sys=0; sys void -FESystem::build_face_table() +FESystem::build_face_tables() { - // Initialize index table - // Multi-component base elements have - // to be thought of. + // Initialize index tables. do this + // in the same way as done for the + // cell tables, except that we now + // loop over the objects of faces + + // For non-primitive shape + // functions, have a special + // invalid index + const std::pair + non_primitive_index (static_cast(-1), + static_cast(-1)); // 1. Vertices unsigned int total_index = 0; - for (unsigned int vertex_number= 0 ; vertex_number < GeometryInfo::vertices_per_face ; + for (unsigned int vertex_number=0; + vertex_number::vertices_per_face; ++vertex_number) { unsigned int comp_start = 0; - for(unsigned int base = 0; base < n_base_elements() ; - ++base) - { - for (unsigned int m = 0; m < element_multiplicity(base); ++m) + for(unsigned int base=0; base - (comp_start+m, - vertex_number*base_element(base).dofs_per_vertex+local_index); + const unsigned int comp_in_base + = base_element(base).face_system_to_component_index(face_index_in_base).first; + const unsigned int comp + = comp_start + comp_in_base; + const unsigned int face_index_in_comp + = base_element(base).face_system_to_component_index(face_index_in_base).second; + face_system_to_component_table[total_index] + = std::make_pair (comp, face_index_in_comp); } + else + face_system_to_component_table[total_index] = non_primitive_index; } - comp_start += element_multiplicity(base); - } } - Assert (total_index <= face_system_to_component_table.size(), - ExcInternalError()); // 2. Lines - for (unsigned line_number= 0 ; ((line_number != GeometryInfo::lines_per_face) && - (GeometryInfo::lines_per_cell > 0)); - ++line_number) - { - unsigned comp_start = 0; - for(unsigned base = 0; base < n_base_elements() ; - ++base) - { - for (unsigned m = 0; m < element_multiplicity(base); ++m) - { - for (unsigned local_index = 0 ; - local_index < base_element(base).dofs_per_line ; - ++local_index) - { - face_system_to_component_table[total_index++] - = std::pair - (comp_start+m, - line_number*base_element(base).dofs_per_line - +local_index+base_element(base).first_face_line_index); - } - } - comp_start += element_multiplicity(base); - } - } - Assert (total_index <= face_system_to_component_table.size(), - ExcInternalError()); + if (GeometryInfo::lines_per_face > 0) + for (unsigned line_number= 0; + line_number != GeometryInfo::lines_per_face; + ++line_number) + { + unsigned comp_start = 0; + for(unsigned base = 0; base < n_base_elements(); ++base) + for (unsigned m=0; m::quads_per_face) && - (GeometryInfo::quads_per_cell > 0)); - ++quad_number) - { - unsigned comp_start = 0; - for(unsigned base = 0; base < n_base_elements() ; - ++base) - { - for (unsigned m = 0; m < element_multiplicity(base); ++m) - { - for (unsigned local_index = 0 ; - local_index < base_element(base).dofs_per_quad ; - ++local_index) - { - face_system_to_component_table[total_index++] - = std::pair - (comp_start+m, - quad_number*base_element(base).dofs_per_quad - +local_index+base_element(base).first_face_quad_index); - } - } - comp_start += element_multiplicity(base); - } - } - Assert (total_index <= face_system_to_component_table.size(), + if (GeometryInfo::quads_per_face > 0) + for (unsigned quad_number= 0; + quad_number != GeometryInfo::quads_per_face; + ++quad_number) + { + unsigned comp_start = 0; + for(unsigned base=0; base dofs_per_component (n_components(), 0); + for (unsigned int sys=0; sys void FESystem::build_interface_constraints () { - // the layout of the constraints matrix is - // described in the FiniteElement class. you - // may want to look there first before trying - // to understand the following, especially - // the mapping of the @p{n} index. +//TODO: Review for base elements with more than one component. do it in the same way as for building the cell restriction and interpolation tables, i.e. top-down + + // the layout of the constraints + // matrix is described in the + // FiniteElement class. you may + // want to look there first before + // trying to understand the + // following, especially the + // mapping of the @p{m} index. // - // in order to map it to the fe-system class, - // we have to know which base element a - // degree of freedom within a vertex, line, - // etc belongs to. this can be accomplished - // by the system_to_component_index + // in order to map it to the + // fe-system class, we have to know + // which base element a degree of + // freedom within a vertex, line, + // etc belongs to. this can be + // accomplished by the + // system_to_component_index // function in conjunction with the - // numbers first_{line,quad,...}_index + // numbers + // first_{line,quad,...}_index for (unsigned int n=0; n n_index - = face_system_to_component_index (n); - - std::pair m_index; + // first for the n + // index. this is simple + // since the n indices are in + // the same order as they are + // usually on a face. note + // that for the data type, + // first value in pair is + // (base element,instance of + // base element), second is + // index within this instance + const std::pair, unsigned int> n_index + = face_system_to_base_table[n]; + + // likewise for the m + // index. this is more + // complicated due to the + // strange ordering we have + // for the dofs on the + // refined faces. + std::pair, unsigned int> m_index; switch (dim) { case 1: @@ -946,12 +1240,12 @@ void FESystem::build_interface_constraints () // as for the first vertex // of the whole cell, so we // can use the - // system_to_component_index - // function (using the - // face_s_t_c_i function would + // system_to_base_table + // variable (using the + // face_s_t_base_t function would // yield the same) if (m < dofs_per_vertex) - m_index = system_to_component_index (m); + m_index = system_to_base_table[m]; else // then come the two sets of // line indices @@ -961,33 +1255,70 @@ void FESystem::build_interface_constraints () const unsigned int sub_line = (m-dofs_per_vertex) / dofs_per_line; Assert (sub_line < 2, ExcInternalError()); - - // get the component by - // asking s_t_c_index and - // tweaking the index a bit - m_index.first = system_to_component_index - (GeometryInfo<2>::vertices_per_cell * dofs_per_vertex - + index_in_line).first; - // first find out the how-many'th - // line index of that component - // this was - m_index.second = (system_to_component_index - (GeometryInfo<2>::vertices_per_cell * dofs_per_vertex - + index_in_line).second - - base_element (component_to_base_table[m_index.first]).first_line_index) - // then add the number of dofs - // per vertex to get the index - // on the first line - + base_element(component_to_base_table[m_index.first]).dofs_per_vertex - // if on the second line: add - // some more - + base_element(component_to_base_table[m_index.first]).dofs_per_line * sub_line; + + // from this + // information, try + // to get base + // element and + // instance of base + // element. we do + // so by + // constructing the + // corresponding + // face index of m + // in the present + // element, then + // use + // face_system_to_base_table + const unsigned int tmp1 = 2*dofs_per_vertex+index_in_line; + m_index.first = face_system_to_base_table[tmp1].first; + + // what we are + // still missing is + // the index of m + // within the base + // elements + // interface_constraints + // table + // + // here, the second + // value of + // face_system_to_base_table + // can help: it + // denotes the face + // index of that + // shape function + // within the base + // element. since + // we know that it + // is a line dof, + // we can construct + // the rest: tmp2 + // will denote the + // index of this + // shape function + // among the line + // shape functions: + Assert (face_system_to_base_table[tmp1].second >= + 2*base_element(m_index.first.first).dofs_per_vertex, + ExcInternalError()); + const unsigned int tmp2 = face_system_to_base_table[tmp1].second - + 2*base_element(m_index.first.first).dofs_per_vertex; + Assert (tmp2 < base_element(m_index.first.first).dofs_per_line, + ExcInternalError()); + m_index.second = base_element(m_index.first.first).dofs_per_vertex + + base_element(m_index.first.first).dofs_per_line*sub_line + + tmp2; }; break; }; case 3: { + // same way as above, + // although a little + // more complicated... + // the indices // m=0..5*d_v-1 are // from the center and @@ -999,76 +1330,65 @@ void FESystem::build_interface_constraints () // so we can use the // simple arithmetic if (m < 5*dofs_per_vertex) - { - m_index.first = system_to_component_index(m % dofs_per_vertex).first; - m_index.second = m / dofs_per_vertex; - } + m_index = system_to_base_table[m]; else // then come the 12 sets of // line indices if (m < 5*dofs_per_vertex + 12*dofs_per_line) - { + { + // for the + // meaning of all + // this, see the + // 2d part const unsigned int index_in_line = (m-5*dofs_per_vertex) % dofs_per_line; const unsigned int sub_line = (m-5*dofs_per_vertex) / dofs_per_line; Assert (sub_line < 12, ExcInternalError()); - - // get the component by - // asking s_t_c_index and - // tweaking the index a bit - m_index.first = system_to_component_index - (GeometryInfo<3>::vertices_per_cell * dofs_per_vertex - + index_in_line).first; - - // first find out the how-many'th - // line index of that component - // this was - m_index.second = (system_to_component_index - (GeometryInfo<3>::vertices_per_cell * dofs_per_vertex - + index_in_line).second - - base_element (component_to_base_table[m_index.first]).first_line_index) - // then add the number of dofs - // for the five vertices to get - // the index on the first line - + 5*base_element(component_to_base_table[m_index.first]).dofs_per_vertex - // and correct for the - // how-many'th line - + base_element(component_to_base_table[m_index.first]).dofs_per_line * sub_line; + + const unsigned int tmp1 = 4*dofs_per_vertex+index_in_line; + m_index.first = face_system_to_base_table[tmp1].first; + + Assert (face_system_to_base_table[tmp1].second >= + 4*base_element(m_index.first.first).dofs_per_vertex, + ExcInternalError()); + const unsigned int tmp2 = face_system_to_base_table[tmp1].second - + 4*base_element(m_index.first.first).dofs_per_vertex; + Assert (tmp2 < base_element(m_index.first.first).dofs_per_line, + ExcInternalError()); + m_index.second = base_element(m_index.first.first).dofs_per_vertex + + base_element(m_index.first.first).dofs_per_line*sub_line + + tmp2; } else // on one of the four sub-quads { + // for the + // meaning of all + // this, see the + // 2d part const unsigned int index_in_quad = (m-5*dofs_per_vertex-12*dofs_per_line) % dofs_per_line; const unsigned int sub_quad = (m-5*dofs_per_vertex-12*dofs_per_line) / dofs_per_line; Assert (sub_quad < 4, ExcInternalError()); - - // get the component by - // asking s_t_c_index and - // tweaking the index a bit - m_index.first = system_to_component_index - (GeometryInfo<3>::vertices_per_cell * dofs_per_vertex - + GeometryInfo<3>::lines_per_cell * dofs_per_line - + index_in_quad).first; - - // first find out the how-many'th - // quad index of that component - // this was - m_index.second = (system_to_component_index - (GeometryInfo<3>::vertices_per_cell * dofs_per_vertex - + GeometryInfo<3>::lines_per_cell * dofs_per_line - + index_in_quad).second - - base_element (component_to_base_table[m_index.first]).first_quad_index) - // then add the number of dofs - // for the five vertices and 12 lines - // to get the index on the first quad - + 5*base_element(component_to_base_table[m_index.first]).dofs_per_vertex - + 12*base_element(component_to_base_table[m_index.first]).dofs_per_line - // and correct for the - // how-many'th line - + base_element(component_to_base_table[m_index.first]).dofs_per_quad * sub_quad; + + const unsigned int tmp1 = 4*dofs_per_vertex+4*dofs_per_line+index_in_quad; + m_index.first = face_system_to_base_table[tmp1].first; + + Assert (face_system_to_base_table[tmp1].second >= + 4*base_element(m_index.first.first).dofs_per_vertex + + 4*base_element(m_index.first.first).dofs_per_line, + ExcInternalError()); + const unsigned int tmp2 = face_system_to_base_table[tmp1].second - + 4*base_element(m_index.first.first).dofs_per_vertex - + 4*base_element(m_index.first.first).dofs_per_line; + Assert (tmp2 < base_element(m_index.first.first).dofs_per_quad, + ExcInternalError()); + m_index.second = base_element(m_index.first.first).dofs_per_vertex + + 12*base_element(m_index.first.first).dofs_per_line + + base_element(m_index.first.first).dofs_per_quad*sub_quad + + tmp2; }; break; @@ -1078,16 +1398,18 @@ void FESystem::build_interface_constraints () Assert (false, ExcNotImplemented()); }; - // now that we gathered all information: - // use it to build the matrix. note - // that if n and m belong to different - // components, there definitely will be - // no coupling + // now that we gathered all + // information: use it to + // build the matrix. note + // that if n and m belong to + // different base elements or + // instances, then there + // definitely will be no + // coupling if (n_index.first == m_index.first) interface_constraints(m,n) - = (base_element(component_to_base_table[n_index.first]) - .constraints()(m_index.second, - n_index.second)); + = (base_element(n_index.first.first).constraints()(m_index.second, + n_index.second)); }; }; @@ -1096,8 +1418,8 @@ void FESystem::build_interface_constraints () template void FESystem::initialize () { - build_cell_table(); - build_face_table(); + build_cell_tables(); + build_face_tables(); // Check if some of the matrices of // the base elements are void. @@ -1122,28 +1444,64 @@ void FESystem::initialize () for (unsigned int i=0;i::children_per_cell;++i) prolongation[i].reinit(0,0); - // distribute the matrices of the base - // finite elements to the matrices of - // this object - for (unsigned int component=0; component::children_per_cell; ++child) { if (do_restriction) - restriction[child] (component_to_system_index (component,i), - component_to_system_index (component, j)) - = base_element(component_to_base_table[component]).restrict(child)(i,j); + restriction[child] (i,j) + = base_element(base).restrict(child)(base_index_i,base_index_j); + if (do_prolongation) - prolongation[child] (component_to_system_index (component,i), - component_to_system_index (component, j)) - = base_element(component_to_base_table[component]).prolongate(child)(i,j); + prolongation[child] (i,j) + = base_element(base).prolongate(child)(base_index_i,base_index_j); }; - + }; // now set up the interface constraints. // this is kind'o hairy, so don't try @@ -1188,7 +1546,7 @@ initialize_unit_support_points () // an empty array to // demonstrate that // fact - for (unsigned int base_el=0 ; base_el > - & base_unit_support_points = base_element(base_el).get_unit_support_points (); - - // otherwise distribute the - // support points of this base - // element to all degrees of - // freedom contributed by this - // base element - Assert(base_unit_support_points.size()==base_element_dofs_per_cell, - ExcInternalError()); - for (unsigned int n=0; n:: initialize_unit_face_support_points () { - // if one of the base elements - // has no support points, then - // it makes no sense to define - // support points for the - // composed element, so return - // an empty array to - // demonstrate that fact (note - // that we ask whether the base - // element has no support - // points at all, not only none - // on the face!) - for (unsigned int base_el=0 ; base_el > & - base_unit_support_points = base_element(base_el).get_unit_face_support_points (); + Assert (index_in_base < base_element(base_i).unit_face_support_points.size(), + ExcInternalError()); - // distribute the support - // points of this base element - // to all degrees of freedom - // contributed by this base - // element - Assert(base_unit_support_points.size()==base_element_dofs_per_face, - ExcNotImplemented()); - for (unsigned int n=0; n::compute_restriction_is_additive_flags (const FiniteElement & +template +std::vector > +FESystem::compute_nonzero_components (const FiniteElement &fe1, + const unsigned int N1) +{ + std::vector*> fe_list; + std::vector multiplicities; + + fe_list.push_back (&fe1); + multiplicities.push_back (N1); + + return compute_nonzero_components (fe_list, multiplicities); +}; + + + +template +std::vector > +FESystem::compute_nonzero_components (const FiniteElement &fe1, + const unsigned int N1, + const FiniteElement &fe2, + const unsigned int N2) +{ + std::vector*> fe_list; + std::vector multiplicities; + + fe_list.push_back (&fe1); + multiplicities.push_back (N1); + + fe_list.push_back (&fe2); + multiplicities.push_back (N2); + + return compute_nonzero_components (fe_list, multiplicities); +}; + + + +template +std::vector > +FESystem::compute_nonzero_components (const FiniteElement &fe1, + const unsigned int N1, + const FiniteElement &fe2, + const unsigned int N2, + const FiniteElement &fe3, + const unsigned int N3) +{ + std::vector*> fe_list; + std::vector multiplicities; + + fe_list.push_back (&fe1); + multiplicities.push_back (N1); + + fe_list.push_back (&fe2); + multiplicities.push_back (N2); + + fe_list.push_back (&fe3); + multiplicities.push_back (N3); + + return compute_nonzero_components (fe_list, multiplicities); +}; + + + +template +std::vector > +FESystem:: +compute_nonzero_components (const std::vector*> &fes, + const std::vector &multiplicities) +{ + Assert (fes.size() == multiplicities.size(), ExcInternalError()); + + // first count the number of dofs + // and components that will emerge + // from the given FEs + unsigned int n_shape_functions = 0; + for (unsigned int i=0; idofs_per_cell * multiplicities[i]; + + unsigned int n_components = 0; + for (unsigned int i=0; in_components() * multiplicities[i]; + + // generate the array that will + // hold the output + std::vector > + retval (n_shape_functions, std::vector (n_components, false)); + + // finally go through all the shape + // functions of the base elements, + // and copy their flags. this + // somehow copies the code in + // build_cell_table, which is not + // nice as it uses too much + // implicit knowledge about the + // layout of the individual bases + // in the composed FE, but there + // seems no way around... + // + // for each shape function, copy + // the non-zero flags from the base + // element to this one, taking into + // account multiplicities, multiple + // components in base elements, and + // other complications + unsigned int total_index = 0; + for (unsigned int vertex_number=0; + vertex_number::vertices_per_cell; + ++vertex_number) + { + unsigned comp_start = 0; + for(unsigned int base=0; basen_components()) + for (unsigned int local_index = 0; + local_index < fes[base]->dofs_per_vertex; + ++local_index, ++total_index) + { + const unsigned int index_in_base + = (fes[base]->dofs_per_vertex*vertex_number + + local_index); + + Assert (comp_start+fes[base]->n_components() <= + retval[total_index].size(), + ExcInternalError()); + for (unsigned int c=0; cn_components(); ++c) + { + Assert (index_in_base < fes[base]->nonzero_components.size(), + ExcInternalError()); + Assert (c < fes[base]->nonzero_components[index_in_base].size(), + ExcInternalError()); + retval[total_index][comp_start+c] + = fes[base]->nonzero_components[index_in_base][c]; + }; + } + } + + // 2. Lines + if (GeometryInfo::lines_per_cell > 0) + for (unsigned int line_number= 0; + line_number != GeometryInfo::lines_per_cell; + ++line_number) + { + unsigned comp_start = 0; + for (unsigned int base=0; basen_components()) + for (unsigned int local_index = 0; + local_index < fes[base]->dofs_per_line; + ++local_index, ++total_index) + { + const unsigned int index_in_base + = (fes[base]->dofs_per_line*line_number + + local_index + + fes[base]->first_line_index); + + Assert (comp_start+fes[base]->n_components() <= + retval[total_index].size(), + ExcInternalError()); + for (unsigned int c=0; cn_components(); ++c) + { + Assert (index_in_base < fes[base]->nonzero_components.size(), + ExcInternalError()); + Assert (c < fes[base]->nonzero_components[index_in_base].size(), + ExcInternalError()); + retval[total_index][comp_start+c] + = fes[base]->nonzero_components[index_in_base][c]; + }; + } + } + + // 3. Quads + if (GeometryInfo::quads_per_cell > 0) + for (unsigned int quad_number= 0; + quad_number != GeometryInfo::quads_per_cell; + ++quad_number) + { + unsigned int comp_start = 0; + for (unsigned int base=0; basen_components()) + for (unsigned int local_index = 0; + local_index < fes[base]->dofs_per_quad; + ++local_index, ++total_index) + { + const unsigned int index_in_base + = (fes[base]->dofs_per_quad*quad_number + + local_index + + fes[base]->first_quad_index); + + Assert (comp_start+fes[base]->n_components() <= + retval[total_index].size(), + ExcInternalError()); + for (unsigned int c=0; cn_components(); ++c) + { + Assert (index_in_base < fes[base]->nonzero_components.size(), + ExcInternalError()); + Assert (c < fes[base]->nonzero_components[index_in_base].size(), + ExcInternalError()); + retval[total_index][comp_start+c] + = fes[base]->nonzero_components[index_in_base][c]; + }; + } + } + + // 4. Hexes + if (GeometryInfo::hexes_per_cell > 0) + for (unsigned int hex_number= 0; + hex_number != GeometryInfo::hexes_per_cell; + ++hex_number) + { + unsigned int comp_start = 0; + for(unsigned int base=0; basen_components()) + for (unsigned int local_index = 0; + local_index < fes[base]->dofs_per_hex; + ++local_index, ++total_index) + { + const unsigned int index_in_base + = (fes[base]->dofs_per_hex*hex_number + + local_index + + fes[base]->first_hex_index); + + Assert (comp_start+fes[base]->n_components() <= + retval[total_index].size(), + ExcInternalError()); + for (unsigned int c=0; cn_components(); ++c) + { + Assert (index_in_base < fes[base]->nonzero_components.size(), + ExcInternalError()); + Assert (c < fes[base]->nonzero_components[index_in_base].size(), + ExcInternalError()); + retval[total_index][comp_start+c] + = fes[base]->nonzero_components[index_in_base][c]; + }; + } + } + + Assert (total_index == retval.size(), ExcInternalError()); + + return retval; +}; + + + template const FiniteElement & FESystem::base_element (const unsigned int index) const @@ -1446,7 +2005,7 @@ FESystem::has_support_on_face (const unsigned int shape_index, { const std::pair component = system_to_component_index(shape_index); - const unsigned int base = component_to_base(component.first); + const unsigned int base = component_to_base(component.first).first; return base_element(base).has_support_on_face(component.second, face_index); } diff --git a/deal.II/deal.II/source/fe/fe_values.cc b/deal.II/deal.II/source/fe/fe_values.cc index a4c2edab9c..aa94880677 100644 --- a/deal.II/deal.II/source/fe/fe_values.cc +++ b/deal.II/deal.II/source/fe/fe_values.cc @@ -34,26 +34,50 @@ using namespace std; template void -FEValuesData::initialize (const unsigned int n_quadrature_points, - const unsigned int n_shapes, - const UpdateFlags flags) +FEValuesData::initialize (const unsigned int n_quadrature_points, + const FiniteElement &fe, + const UpdateFlags flags) { update_flags = flags; - + + // initialize the table mapping + // from shape function number to + // the rows in the tables denoting + // its first non-zero + // component. with this also count + // the total number of non-zero + // components accumulated over all + // shape functions + shape_function_to_row_table.resize (fe.dofs_per_cell); + unsigned int row = 0; + for (unsigned int i=0; i= fe.dofs_per_cell, + ExcInternalError()); + + // with the number of rows now + // known, initialize those fields + // that we will need to their + // correct size if (flags & update_values) - shape_values.reinit(n_shapes, n_quadrature_points); + shape_values.reinit(n_nonzero_shape_components, n_quadrature_points); if (flags & update_gradients) { - shape_gradients.resize(n_shapes); - for (unsigned int i=0;i::FEValuesBase (const unsigned int n_q_points, template FEValuesBase::~FEValuesBase () { - if (fe_data) + // delete those fields that were + // created by the mapping and + // finite element objects, + // respectively, but of which we + // have assumed ownership + if (fe_data != 0) { typename Mapping::InternalDataBase *tmp1=fe_data; fe_data=0; delete tmp1; } - if (mapping_data) + if (mapping_data != 0) { typename Mapping::InternalDataBase *tmp1=mapping_data; mapping_data=0; @@ -138,11 +167,15 @@ void FEValuesBase::get_function_values (const InputVector &fe_fu std::fill_n (values.begin(), n_quadrature_points, 0); // add up contributions of trial + // functions. note that here we + // deal with scalar finite + // elements, so no need to check + // for non-primitivity of shape // functions for (unsigned int point=0; point::get_function_values (const InputVector { Assert (n_quadrature_points == values.size(), ExcWrongVectorSize(values.size(), n_quadrature_points)); + + const unsigned int n_components = fe->n_components(); for (unsigned i=0;in_components(), - ExcWrongNoOfComponents()); + Assert (values[i].size() == n_components, ExcWrongNoOfComponents()); + Assert (update_flags & update_values, ExcAccessToUninitializedField()); Assert (fe_function.size() == present_cell->get_dof_handler().n_dofs(), ExcWrongVectorSize(fe_function.size(), present_cell->get_dof_handler().n_dofs())); @@ -174,42 +209,20 @@ void FEValuesBase::get_function_values (const InputVector std::fill_n (values[i].begin(), values[i].size(), 0); // add up contributions of trial - // functions + // functions. now check whether the + // shape function is primitive or + // not. if it is, then set its only + // non-zero component, otherwise + // loop over components for (unsigned int point=0; pointsystem_to_component_index(shape_func).first) - += (dof_values(shape_func) * shape_values(shape_func, point)); -}; - - - -template -const typename FEValuesData::ShapeVector & -FEValuesBase::get_shape_values () const -{ - Assert (update_flags & update_values, ExcAccessToUninitializedField()); - return shape_values; -}; - - - -template -const typename FEValuesData::GradientVector & -FEValuesBase::get_shape_grads () const -{ - Assert (update_flags & update_gradients, ExcAccessToUninitializedField()); - return shape_gradients; -}; - - - - -template -const typename std::vector > > & -FEValuesBase::get_shape_2nd_derivatives () const -{ - Assert (update_flags & update_second_derivatives, ExcAccessToUninitializedField()); - return shape_2nd_derivatives; + if (fe->is_primitive(shape_func)) + values[point](fe->system_to_component_index(shape_func).first) + += (dof_values(shape_func) * shape_value(shape_func, point)); + else + for (unsigned int c=0; c()); // add up contributions of trial + // functions. note that here we + // deal with scalar finite + // elements, so no need to check + // for non-primitivity of shape // functions for (unsigned int point=0; point tmp(shape_gradients[shape_func][point]); + Tensor<1,dim> tmp = shape_grad(shape_func,point); tmp *= dof_values(shape_func); gradients[point] += tmp; }; @@ -283,9 +300,11 @@ get_function_grads (const InputVector { Assert (n_quadrature_points == gradients.size(), ExcWrongNoOfComponents()); - for (unsigned i=0;in_components(), - ExcWrongVectorSize(gradients[i].size(), fe->n_components())); + + const unsigned int n_components = fe->n_components(); + for (unsigned i=0; iget_dof_handler().n_dofs(), ExcWrongVectorSize(fe_function.size(), present_cell->get_dof_handler().n_dofs())); @@ -303,15 +322,27 @@ get_function_grads (const InputVector std::fill_n (gradients[i].begin(), gradients[i].size(), Tensor<1,dim>()); // add up contributions of trial - // functions + // functions. now check whether the + // shape function is primitive or + // not. if it is, then set its only + // non-zero component, otherwise + // loop over components for (unsigned int point=0; point tmp(shape_gradients[shape_func][point]); - tmp *= dof_values(shape_func); - gradients[point][fe->system_to_component_index(shape_func).first] - += tmp; - }; + if (fe->is_primitive (shape_func)) + { + Tensor<1,dim> tmp = shape_grad(shape_func,point); + tmp *= dof_values(shape_func); + gradients[point][fe->system_to_component_index(shape_func).first] + += tmp; + } + else + for (unsigned int c=0; c tmp = shape_grad_component(shape_func,point,c); + tmp *= dof_values(shape_func); + gradients[point][c] += tmp; + }; }; @@ -343,11 +374,15 @@ get_function_2nd_derivatives (const InputVector &fe_function, std::fill_n (second_derivatives.begin(), n_quadrature_points, Tensor<2,dim>()); // add up contributions of trial + // functions. note that here we + // deal with scalar finite + // elements, so no need to check + // for non-primitivity of shape // functions for (unsigned int point=0; point tmp(shape_2nd_derivatives[shape_func][point]); + Tensor<2,dim> tmp = shape_2nd_derivative(shape_func,point); tmp *= dof_values(shape_func); second_derivatives[point] += tmp; }; @@ -364,9 +399,11 @@ get_function_2nd_derivatives (const InputVector { Assert (n_quadrature_points == second_derivs.size(), ExcWrongNoOfComponents()); + + const unsigned int n_components = fe->n_components(); for (unsigned i=0;in_components(), - ExcWrongVectorSize(second_derivs[i].size(), fe->n_components())); + Assert (second_derivs[i].size() == n_components, ExcWrongNoOfComponents()); + Assert (update_flags & update_second_derivatives, ExcAccessToUninitializedField()); Assert (fe_function.size() == present_cell->get_dof_handler().n_dofs(), ExcWrongVectorSize(fe_function.size(), present_cell->get_dof_handler().n_dofs())); @@ -387,12 +424,20 @@ get_function_2nd_derivatives (const InputVector // functions for (unsigned int point=0; point tmp(shape_2nd_derivatives[shape_func][point]); - tmp *= dof_values(shape_func); - second_derivs[point][fe->system_to_component_index(shape_func).first] - += tmp; - }; + if (fe->is_primitive(shape_func)) + { + Tensor<2,dim> tmp(shape_2nd_derivative(shape_func,point)); + tmp *= dof_values(shape_func); + second_derivs[point][fe->system_to_component_index(shape_func).first] + += tmp; + } + else + for (unsigned int c=0; c tmp = shape_2nd_derivative_component(shape_func,point,c); + tmp *= dof_values(shape_func); + second_derivs[point][c] += tmp; + }; }; @@ -441,7 +486,8 @@ FEValuesBase::memory_consumption () const MemoryConsumption::memory_consumption (mapping_data) + MemoryConsumption::memory_consumption (*mapping_data) + MemoryConsumption::memory_consumption (fe_data) + - MemoryConsumption::memory_consumption (*fe_data)); + MemoryConsumption::memory_consumption (*fe_data) + + MemoryConsumption::memory_consumption (shape_function_to_row_table)); }; @@ -536,9 +582,7 @@ FEValues::initialize (const UpdateFlags update_flags) fe_data = fe->get_data(flags, *mapping, quadrature); // set up objects within this class - FEValuesData::initialize(n_quadrature_points, - dofs_per_cell, - flags); + FEValuesData::initialize (n_quadrature_points, *fe, flags); }; @@ -685,9 +729,7 @@ FEFaceValues::initialize (const UpdateFlags update_flags) fe_data = fe->get_face_data(flags, *mapping, quadrature); // set up objects within this class - FEValuesData::initialize(n_quadrature_points, - dofs_per_cell, - flags); + FEValuesData::initialize(n_quadrature_points, *fe, flags); }; @@ -775,9 +817,7 @@ FESubfaceValues::initialize (const UpdateFlags update_flags) fe_data = fe->get_subface_data(flags, *mapping, quadrature); // set up objects within this class - FEValuesData::initialize(n_quadrature_points, - dofs_per_cell, - flags); + FEValuesData::initialize(n_quadrature_points, *fe, flags); }; diff --git a/deal.II/doc/news/2002/c-3-4.html b/deal.II/doc/news/2002/c-3-4.html index 4ca67185d3..f3513fcdaa 100644 --- a/deal.II/doc/news/2002/c-3-4.html +++ b/deal.II/doc/news/2002/c-3-4.html @@ -33,15 +33,47 @@ contributor's names are abbreviated by WB (Wolfgang Bangerth), GK

General

    +
  1. + Removed: The functions FEValuesBase::get_shape_values, FEValuesBase::get_shape_grads, and FEValuesBase::get_shape_2nd_derivatives are + now removed as they expose too much of the internal data + structure of their class, and interfere with plans for the + extension of this and related classes. These functions, which + had been deprecated in the previous version, are thus not + part of the release after version 3.3 any more. Use the FEValuesBase::shape_value and alike + functions as a replacement. +
    + For more information, read + this mail. +
    + (WB 2002/06/10) +

    -
  2. Improved: make knows a target veryclean - now. clean leaves the libraries now, removing everything - that is not needed to use - deal.II. veryclean removes even the - libraries, leaving the directory more or less in the state like - after configure. -
    - (GK 2002/06/07) +
  3. + New: deal.II now also supports vector-valued + finite elements with shape functions for which more than just + one vector component is non-zero. Such elements are, for + example, the Nedelec and Raviart-Thomas families. Previously, + vector-valued elements were only supported insofar as they + could be composed of scalar ones; in that case, each + (vector-valued) shape function had only one non-zero vector + component. +
    + (WB 2002/06/10) +

    + +
  4. + Improved: make knows a target veryclean + now. clean leaves the libraries now, removing + everything that is not needed to use + deal.II. veryclean removes even the + libraries, leaving the directory more or less in the state like + after configure. +
    + (GK 2002/06/07)

diff --git a/deal.II/doc/news/news.html b/deal.II/doc/news/news.html index a64502a5d9..5ee6713996 100644 --- a/deal.II/doc/news/news.html +++ b/deal.II/doc/news/news.html @@ -48,6 +48,21 @@
+
+ 2002/06/10: Supporting framework for vector-valued + finite elements committed +
+
+ deal.II supports vector-valued finite + elements for a long time, but only by composing them of + scalar elements. This excludes Nedelec or Raviart-Thomas + type elements, for which every shape function is + vector-valued and non-zero in all components. Support for + such elements is now available, and will be part of the next + version. +

+
+
Changes between version 3.4 and the main branch diff --git a/tests/fe/Makefile b/tests/fe/Makefile index 1abc522555..c4470a7eb9 100644 --- a/tests/fe/Makefile +++ b/tests/fe/Makefile @@ -24,15 +24,18 @@ default: run-tests derivatives.exe : derivatives.go $(libraries) fe_data_test.exe : fe_data_test.go $(libraries) +fe_traits.exe : fe_traits.go $(libraries) +internals.exe : internals.go $(libraries) mapping.exe : mapping.go $(libraries) mapping_c1.exe : mapping_c1.go $(libraries) -shapes.exe : shapes.go $(libraries) -numbering.exe : numbering.go $(libraries) mapping_q1_eulerian.exe : mapping_q1_eulerian.go $(libraries) -internals.exe : internals.go $(libraries) +non_primitive_1.exe : non_primitive_1.go $(libraries) +numbering.exe : numbering.go $(libraries) +shapes.exe : shapes.go $(libraries) -tests = derivatives fe_data_test mapping mapping_c1 shapes \ - numbering mapping_q1_eulerian internals +tests = derivatives fe_data_test fe_traits internals \ + mapping mapping_c1 mapping_q1_eulerian \ + non_primitive_1 numbering shapes ############################################################ diff --git a/tests/fe/fe_data_test.cc b/tests/fe/fe_data_test.cc index 73a082b77a..41c28cc7b8 100644 --- a/tests/fe/fe_data_test.cc +++ b/tests/fe/fe_data_test.cc @@ -2,7 +2,7 @@ // $Id$ // Version: $Name$ // -// Copyright (C) 1998, 1999, 2000, 2001 by the deal.II authors +// Copyright (C) 1998, 1999, 2000, 2001, 2002 by the deal.II authors // // This file is subject to QPL and may not be distributed // without copyright and license information. Please refer @@ -19,6 +19,7 @@ #include #include #include +#include #include @@ -40,6 +41,22 @@ void test_fe_datas() FE_Q (3), 1)); fe_datas.push_back(new FESystem(FE_Q (4), 2)); } + + // have systems of systems, and + // construct hierarchies of + // subsequently weirder elements by + // taking each of them in turn as + // basis of others + fe_datas.push_back (new FESystem (FESystem (FE_Q(1),2),2)); + fe_datas.push_back (new FESystem (FESystem (FE_Q(1),2),1, + FESystem (FE_DGQ(1),2),1)); + fe_datas.push_back (new FESystem (FESystem (FE_Q(1),1, + FE_Q(2),1),1, + FESystem (FE_Q(2),2),1, + FESystem (FE_DGQ(2),2),1)); + fe_datas.push_back (new FESystem (*fe_datas[fe_datas.size()-3], 2, + *fe_datas[fe_datas.size()-2], 1, + *fe_datas[fe_datas.size()-1], 2)); deallog << std::endl << "dim=" << dim << std::endl; for (unsigned int n=0; n +#include + +#include +#include +#include +#include +#include +#include +#include + + +template +void check (const FiniteElement &fe) +{ + const unsigned int dofs_per_cell = fe.dofs_per_cell; + + // first check whether shape + // functions are primitive: + deallog << " Primitivity: "; + for (unsigned int i=0; i +void check () +{ + // check usual Lagrange elements + for (unsigned int p=1; p<3; ++p) + { + deallog << "Checking FE_Q<" << dim << ">(" << p << "): " + << std::endl; + check (FE_Q(p)); + }; + + // check DG Lagrange elements + for (unsigned int p=0; p<3; ++p) + { + deallog << "Checking FE_DGQ<" << dim << ">(" << p << "): " + << std::endl; + check (FE_DGQ(p)); + }; + + // check DG-P elements + for (unsigned int p=0; p<3; ++p) + { + deallog << "Checking FE_DGP<" << dim << ">(" << p << "): " + << std::endl; + check (FE_DGP(p)); + }; + + // check systems of Q-elements + for (unsigned int p=1; p<3; ++p) + { + deallog << "Checking FE_Q<" << dim << ">(" << p << ")^2: " + << std::endl; + check (FESystem (FE_Q(p),2)); + }; + + // check systems of systems of + // Q-elements + for (unsigned int p=1; p<3; ++p) + { + deallog << "Checking FE_Q<" << dim << ">(" << p << ")^2^2: " + << std::endl; + check (FESystem (FESystem (FE_Q(p),2), 2)); + }; +}; + + + + +int main () +{ + std::ofstream logfile("fe_traits.output"); + deallog.attach(logfile); + deallog.depth_console(0); + + check<1> (); + check<2> (); + check<3> (); +}; + diff --git a/tests/fe/shapes.cc b/tests/fe/shapes.cc index 3956cfbf7e..f8ccf43bcc 100644 --- a/tests/fe/shapes.cc +++ b/tests/fe/shapes.cc @@ -39,7 +39,9 @@ plot_shape_functions(Mapping& mapping, QTrapez<1> q_trapez; QIterated q(q_trapez, div); - FEValues fe(mapping, finel, q, UpdateFlags(update_values)); + FEValues fe(mapping, finel, q, UpdateFlags(update_values | + update_gradients | + update_second_derivatives)); fe.reinit(c); sprintf(fname, "Shapes%dd-%s", dim, name); @@ -57,6 +59,27 @@ plot_shape_functions(Mapping& mapping, for (unsigned int i=0;i()) && + (fe.shape_2nd_derivative_component(i,k,c) == Tensor<2,dim>()), + ExcInternalError()); + }; } deallog << std::endl; k++; @@ -93,9 +116,13 @@ plot_face_shape_functions(Mapping& mapping, QTrapez<1> q_trapez; QIterated q(q_trapez, div); FEFaceValues fe(mapping, finel, q, UpdateFlags(update_values - | update_q_points)); + | update_gradients + | update_second_derivatives + | update_q_points)); FESubfaceValues sub(mapping, finel, q, UpdateFlags(update_values - | update_q_points)); + | update_gradients + | update_second_derivatives + | update_q_points)); sprintf(fname, "ShapesFace%dd-%s", dim, name); deallog.push(fname); @@ -116,6 +143,27 @@ plot_face_shape_functions(Mapping& mapping, for (unsigned int i=0;i()) && + (fe.shape_2nd_derivative_component(i,k,c) == Tensor<2,dim>()), + ExcInternalError()); + }; } deallog << std::endl; k++; @@ -138,6 +186,27 @@ plot_face_shape_functions(Mapping& mapping, for (unsigned int i=0;i()) && + (sub.shape_2nd_derivative_component(i,k,c) == Tensor<2,dim>()), + ExcInternalError()); + }; } deallog << std::endl; k++; @@ -178,9 +247,19 @@ void test_compute_functions (const Mapping &mapping, bool coincide=true; for (unsigned int x=0; x1e-14) - coincide=false; + { + if (fabs(fe_values.shape_value(i,x)-fe.shape_value(i,q.point(x)))>1e-14) + coincide=false; + for (unsigned int c=0; c &mapping, tmp-=fe.shape_grad(i,q.point(x)); if (sqrt(tmp*tmp)>1e-14) coincide=false; + + for (unsigned int c=0; c())), + ExcInternalError()); } if (!coincide) @@ -211,8 +298,16 @@ void test_compute_functions (const Mapping &mapping, if (diff>max_diff) max_diff=diff; if (fabs(tmp[j][k])>1e-6) coincide=false; - } - + } + + for (unsigned int c=0; c())), + ExcInternalError()); } if (!coincide)