From: Wolfgang Bangerth Date: Fri, 24 Apr 1998 14:13:33 +0000 (+0000) Subject: Add more support for curved boundaries and the integration on faces. X-Git-Tag: v8.0.0~23063 X-Git-Url: https://gitweb.dealii.org/cgi-bin/gitweb.cgi?a=commitdiff_plain;h=b13a6c153b667e2fbf0530811b58ef6181402084;p=dealii.git Add more support for curved boundaries and the integration on faces. git-svn-id: https://svn.dealii.org/trunk@191 0785d39b-7218-0410-832d-ea1e28bc413d --- diff --git a/deal.II/deal.II/Todo b/deal.II/deal.II/Todo index 859ab7724e..a235c34cee 100644 --- a/deal.II/deal.II/Todo +++ b/deal.II/deal.II/Todo @@ -20,12 +20,17 @@ Write monitors to control whether enough memory was allocated for Write Triangulation::mesh_function. + Unify CellAccessor<1> and <2> by renaming LineAccessor -> SubstructAccessor QuadAccessor -> SubstructAccessor and deriving CellAccessor from SubstructAccessor Do the same with DoFLineAccessor and DoFQuadAccessor +Unify FEValues and FEFaceValues to FEValues + and use partial specialization. + + Add support for fast assemblage of mass matrices. Make AssemblerData a local class to Assembler again if gcc2.8 supports it. diff --git a/deal.II/deal.II/include/fe/fe.h b/deal.II/deal.II/include/fe/fe.h index 482e694caa..8a81163cbb 100644 --- a/deal.II/deal.II/include/fe/fe.h +++ b/deal.II/deal.II/include/fe/fe.h @@ -12,6 +12,7 @@ // forward declarations +template class Boundary; template class FiniteElement; template class Quadrature; @@ -127,6 +128,13 @@ class FEValues { * specified finite element using the * quadrature points of the given * quadrature rule. + * + * This function actually only fills + * the fields related to the unit face, + * the fields related to a real face (like + * gradients, true quadrature points, etc.) + * need to be initialized using the + * #reinit# function. */ FEValues (const FiniteElement &, const Quadrature &, @@ -228,9 +236,23 @@ class FEValues { * Reinitialize the gradients, Jacobi * determinants, etc for the given cell * and the given finite element. + * + * This function needs a boundary object + * passed, since this class needs to know + * how to handle faces which are located + * on the boundary of the domain. In that + * case, faces may be curved and the + * calculation of quadrature points, + * gradients and the like may need + * additional effort, depending on the + * mapping from the unit to the real cell + * (linear mappings use straight boundary + * segments, but higher order elements + * may use other ways.) */ void reinit (const Triangulation::cell_iterator &, - const FiniteElement &); + const FiniteElement &, + const Boundary &); /** * Exception @@ -346,6 +368,306 @@ class FEValues { +/** + Represent a finite element evaluated with a specific quadrature rule. + This class is an optimization which avoids evaluating the shape functions + at the quadrature points each time a quadrature takes place. Rather, the + values and gradients (and possibly higher order derivatives in future + versions of this library) are evaluated once and for all on the unit + face before doing the quadrature itself. Only the Jacobian matrix of + the transformation from the unit face to the real face and the integration + points in real space are calculated each time we move on to a new face. + + The unit face is defined to be the tensor product of the interval $[0,1]$ + in the present number of dimensions minus one. In part of the literature, + the convention is used that the unit cell be the tensor product of the + interval $[-1,1]$, which is to distinguished properly. + + This class is very similar to the #FEValues# class; see there for more + documentation. + */ +template +class FEFaceValues { + public: + + + + /** + * Number of quadrature points on + * the face. + */ + const unsigned int n_quadrature_points; + + /** + * Total number of shape functions + * on the cell adjacent to this face. + * This number is not the same as the + * number of shape functions of which + * the center is located on the face. + */ + const unsigned int total_dofs; + + /** + * Constructor. Fill all arrays with the + * values of the shape functions of the + * specified finite element using the + * quadrature points of the given + * quadrature rule for the face, which + * has a dimension one less than the + * cell. + * + * This function actually only fills + * the fields related to the unit face, + * the fields related to a real face (like + * gradients, true quadrature points, etc.) + * need to be initialized using the + * #reinit# function. + */ + FEFaceValues (const FiniteElement &, + const Quadrature &, + const UpdateFields); + + /** + * Return the value of the #i#th shape + * function at the #j# quadrature point. + */ + double shape_value (const unsigned int i, + const unsigned int j) const; + + /** + * Return a pointer to the matrix holding + * all values 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 dFMatrix & get_shape_values () const; + + /** + * Return the gradient of the #i#th shape + * function at the #j# quadrature point. + * If you want to get the derivative in + * one of the coordinate directions, use + * the appropriate function of the #Point# + * class to extract one component. Since + * only a reference to the gradient's value + * is returned, there should be no major + * performance drawback. + * The function returns the gradient on the + * real element, not the reference element. + */ + const Point & shape_grad (const unsigned int i, + const unsigned int j) const; + + /** + * Return a 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 vector > > & get_shape_grads () const; + + /** + * Return the position of the #i#th + * quadrature point in real space. + * + * For curved boundary cells, using + * biquadratic or higher mappings + * of the unit cell to the real cell, + * these points may not be on the + * plane submannifold on which the + * vertices of the face lie. + */ + const Point & quadrature_point (const unsigned int i) const; + + /** + * Return a pointer to the vector of + * quadrature points. + */ + const vector > & get_quadrature_points () const; + + /** + * Return the point in real space where + * the #i#th ansatz function is located + * (location is in the sense of where it + * assumes its nominal properties, e.g. at + * the vertex of a cell, at the center of + * a line, etc). + * + * This function is needed for the + * interpolation problem: if we want to + * transfer a continuous function to a + * finite element function by interpolation + * we have to take the continuous + * function's value at the ansatz function + * locations. + */ + const Point & ansatz_point (const unsigned int i) const; + + /** + * Return a pointer to the vector of points + * denoting the location of the ansatz + * functions. + */ + const vector > & get_ansatz_points () const; + + /** + * Return the Jacobi determinant times + * the weight of the #i#th quadrature + * point. + */ + double JxW (const unsigned int i) const; + + /** + * Return a pointer to the array holding + * the JxW values at the different + * quadrature points. + */ + const vector & get_JxW_values () const; + + /** + * Reinitialize the gradients, Jacobi + * determinants, etc for the given cell + * and the given finite element. + * + * The constructor needs a boundary object + * passed, since this class needs to know + * how to handle faces which are located + * on the boundary of the domain. In that + * case, faces may be curved and the + * calculation of quadrature points, + * gradients and the like may need + * additional effort, depending on the + * mapping from the unit to the real cell + * (linear mappings use straight boundary + * segments, but higher order elements + * may use other ways.) + */ + void reinit (const Triangulation::face_iterator &, + const FiniteElement &, + const Boundary &); + + /** + * Exception + */ + DeclException2 (ExcInvalidIndex, + int, int, + << "The index " << arg1 + << " is out of range, it should be less than " << arg2); + /** + * Exception + */ + DeclException0 (ExcAccessToUninitializedField); + /** + * Exception + */ + DeclException0 (ExcCannotInitializeField); + + private: + /** + * Store the values of the shape functions + * at the quadrature points. Rows in this + * matrix denote the values of a single + * shape function at the different points, + * columns are for a single point with the + * different shape functions. + */ + dFMatrix shape_values; + + /** + * Store the gradients of the shape + * functions at the quadrature points. + * Since unfortunately the full matrix + * classes of DEAL are not templated, + * we have to store them in an + * archetypic style. + * + * This field is reset each time + * #reinit# is called and contains the + * gradients on the real element, rather + * than on the reference element. + */ + vector > > shape_gradients; + + /** + * Store the gradients of the shape + * functions at the quadrature points on + * the unit cell. + * This field is set up upon construction + * of the object and contains the gradients + * on the reference element. + */ + vector > > unit_shape_gradients; + + /** + * Store an array of the weights of the + * quadrature points. This array is + * set up upon construction. + */ + vector weights; + + /** + * Store an array of weights times the + * Jacobi determinant at the quadrature + * points. This function is reset each time + * #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. + */ + vector JxW_values; + + /** + * Array of quadrature points. This array + * is set up upon calling #reinit# and + * contains the quadrature points on the + * real element, rather than on the + * reference element. + */ + vector > quadrature_points; + + /** + * Array of quadrature points in the unit + * cell. This array is set up upon + * construction and contains the quadrature + * points on the reference element. + */ + vector > unit_quadrature_points; + + /** + * Array of points denoting the off-point + * of the ansatz functions. In real space + * (no-one seems to need the off-point + * on the unit cell, so no function is + * provided for this). + */ + vector > ansatz_points; + + /** + * Store the jacobi matrices at the + * different quadrature points. This field + * is set each time #reinit# is called. + */ + vector jacobi_matrices; + + /** + * Store a pointer to the object describing + * the boundary of the domain. + */ + const Boundary &boundary; + + /** + * Store which fields are to be updated by + * the reinit function. + */ + UpdateFields update_flags; +}; + + + + /** Base class for finite elements in arbitrary dimensions. */ @@ -490,7 +812,8 @@ class FiniteElementBase { vector > &ansatz_points, const bool compute_ansatz_points, vector > &q_points, - const bool compute_q_points) const; + const bool compute_q_points, + const Boundary &boundary) const; /** * Return the ansatz points this FE has @@ -498,13 +821,15 @@ class FiniteElementBase { * given face as a side. This function is * needed for higher order elements, if * we want to use curved boundary - * approximations. + * approximations. For that reason, a + * boundary object has to be passed. * * The function assumes that the fields * already have the right number of * elements. */ virtual void face_ansatz_points (const Triangulation::face_iterator &face, + const Boundary &boundary, vector > &ansatz_points) const; /** @@ -719,7 +1044,8 @@ class FiniteElement<1> : public FiniteElementBase<1> { vector > &ansatz_points, const bool compute_ansatz_points, vector > &q_points, - const bool compute_q_points) const; + const bool compute_q_points, + const Boundary<1> &boundary) const; /** * Return the ansatz points this FE has @@ -733,6 +1059,7 @@ class FiniteElement<1> : public FiniteElementBase<1> { * At present it is not implemented. */ virtual void face_ansatz_points (const Triangulation<1>::face_iterator &face, + const Boundary<1> &boundary, vector > &ansatz_points) const; }; @@ -869,7 +1196,8 @@ class FiniteElement<2> : public FiniteElementBase<2> { vector > &ansatz_points, const bool compute_ansatz_points, vector > &q_points, - const bool compute_q_points) const; + const bool compute_q_points, + const Boundary<2> &boundary) const; /** * Return the ansatz points this FE has @@ -884,6 +1212,7 @@ class FiniteElement<2> : public FiniteElementBase<2> { * elements. */ virtual void face_ansatz_points (const Triangulation<2>::face_iterator &face, + const Boundary<2> &boundary, vector > &ansatz_points) const; }; diff --git a/deal.II/deal.II/include/fe/fe_lib.lagrange.h b/deal.II/deal.II/include/fe/fe_lib.lagrange.h index e4da5ae11c..792afc5073 100644 --- a/deal.II/deal.II/include/fe/fe_lib.lagrange.h +++ b/deal.II/deal.II/include/fe/fe_lib.lagrange.h @@ -80,7 +80,8 @@ class FELinear : public FiniteElement { vector > &ansatz_points, const bool compute_ansatz_points, vector > &q_points, - const bool compute_q_points) const; + const bool compute_q_points, + const Boundary &boundary) const; /** * Return the ansatz points this FE has @@ -95,6 +96,7 @@ class FELinear : public FiniteElement { * elements. */ virtual void face_ansatz_points (const Triangulation::face_iterator &face, + const Boundary &boundary, vector > &ansatz_points) const; }; @@ -164,7 +166,8 @@ class FEQuadratic : public FiniteElement { vector > &ansatz_points, const bool compute_ansatz_points, vector > &q_points, - const bool compute_q_points) const; + const bool compute_q_points, + const Boundary &boundary) const; }; @@ -233,7 +236,8 @@ class FECubic : public FiniteElement { vector > &ansatz_points, const bool compute_ansatz_points, vector > &q_points, - const bool compute_q_points) const; + const bool compute_q_points, + const Boundary &boundary) const; }; diff --git a/deal.II/deal.II/include/numerics/assembler.h b/deal.II/deal.II/include/numerics/assembler.h index f8c5f78888..fa1ad5121f 100644 --- a/deal.II/deal.II/include/numerics/assembler.h +++ b/deal.II/deal.II/include/numerics/assembler.h @@ -106,7 +106,8 @@ struct AssemblerData { dVector &rhs_vector, const Quadrature &quadrature, const FiniteElement &fe, - const UpdateFields &update_flags); + const UpdateFields &update_flags, + const Boundary &boundary); /** * Pointer to the dof handler object @@ -151,7 +152,16 @@ struct AssemblerData { * FEValues object need to be reinitialized * on each cell. */ - const UpdateFields update_flags; + const UpdateFields update_flags; + + /** + * Store a pointer to the object describing + * the boundary of the domain. This is + * necessary, since we may want to use + * curved faces of cells at the boundary + * when using higher order elements. + */ + const Boundary &boundary; }; @@ -229,13 +239,13 @@ class Assembler : public DoFCellAccessor { * Pointer to the matrix to be assembled * by this object. */ - dSMatrix &matrix; + dSMatrix &matrix; /** * Pointer to the vector to be assembled * by this object. */ - dVector &rhs_vector; + dVector &rhs_vector; /** * Pointer to the finite element used for @@ -258,6 +268,15 @@ class Assembler : public DoFCellAccessor { * quadrature points. */ FEValues fe_values; + + /** + * Store a pointer to the object describing + * the boundary of the domain. This is + * necessary, since we may want to use + * curved faces of cells at the boundary + * when using higher order elements. + */ + const Boundary &boundary; }; diff --git a/deal.II/deal.II/include/numerics/base.h b/deal.II/deal.II/include/numerics/base.h index 0b943cfd88..b4cf45b036 100644 --- a/deal.II/deal.II/include/numerics/base.h +++ b/deal.II/deal.II/include/numerics/base.h @@ -20,7 +20,8 @@ template class DataOut; template class Function; template class Equation; template class Assembler; - +template class Boundary; +template class StraightBoundary; @@ -78,6 +79,13 @@ enum NormType { induced by hanging nodes. \end{itemize} + The #assemble# function needs an object describing the boundary of the domain, + since for higher order finite elements, we may be tempted to use curved faces + of cells for better approximation of the boundary. In this case, the + transformation from the unit cell to the real cell requires knowledge of + the exact boundary of the domain. By default it is assumed that the + boundary be approximated by straight segments. + {\bf Solving} @@ -294,7 +302,8 @@ class ProblemBase { const Quadrature &q, const FiniteElement &fe, const UpdateFields &update_flags, - const DirichletBC &dirichlet_bc = DirichletBC()); + const DirichletBC &dirichlet_bc = DirichletBC(), + const Boundary &boundary = StraightBoundary()); /** * Solve the system of equations. @@ -315,7 +324,8 @@ class ProblemBase { dVector &difference, const Quadrature &q, const FiniteElement &fe, - const NormType &norm) const; + const NormType &norm, + const Boundary &boundary=StraightBoundary()) const; /** * Initialize the #DataOut# object with @@ -358,6 +368,7 @@ class ProblemBase { */ virtual void make_boundary_value_list (const DirichletBC &dirichlet_bc, const FiniteElement &fe, + const Boundary &boundary, map &boundary_values) const; /** @@ -433,7 +444,8 @@ class ProblemBase { dVector &solution, dVector &right_hand_side, const DirichletBC &dirichlet_bc, - const FiniteElement &fe); + const FiniteElement &fe, + const Boundary &boundary); friend class Assembler; }; diff --git a/deal.II/deal.II/source/fe/fe.cc b/deal.II/deal.II/source/fe/fe.cc index b28f842584..911f77d438 100644 --- a/deal.II/deal.II/source/fe/fe.cc +++ b/deal.II/deal.II/source/fe/fe.cc @@ -4,7 +4,7 @@ #include #include #include - +#include @@ -52,8 +52,8 @@ FEValues::FEValues (const FiniteElement &fe, template double FEValues::shape_value (const unsigned int i, const unsigned int j) const { - Assert (i<(unsigned int)shape_values.m(), ExcInvalidIndex (i, shape_values.m())); - Assert (j<(unsigned int)shape_values.n(), ExcInvalidIndex (j, shape_values.n())); + Assert (i const Point & FEValues::shape_grad (const unsigned int i, const unsigned int j) const { - Assert (i<(unsigned int)shape_values.m(), ExcInvalidIndex (i, shape_values.m())); - Assert (j<(unsigned int)shape_values.n(), ExcInvalidIndex (j, shape_values.n())); + Assert (i::JxW (const unsigned int i) const { template void FEValues::reinit (const typename Triangulation::cell_iterator &cell, - const FiniteElement &fe) { + const FiniteElement &fe, + const Boundary &boundary) { // fill jacobi matrices and real // quadrature points if ((update_flags | update_jacobians) || @@ -117,7 +118,8 @@ void FEValues::reinit (const typename Triangulation::cell_iterator &ce ansatz_points, update_flags | update_ansatz_points, quadrature_points, - update_flags | update_q_points); + update_flags | update_q_points, + boundary); // compute gradients on real element if // requested @@ -149,8 +151,7 @@ void FEValues::reinit (const typename Triangulation::cell_iterator &ce // determinant if (update_flags | update_JxW_values) { - Assert (update_flags | update_jacobians, - ExcCannotInitializeField()); + Assert (update_flags | update_jacobians, ExcCannotInitializeField()); for (unsigned int i=0; i::fill_fe_values (const typename Triangulation:: vector > &, const bool, vector > &, - const bool) const { + const bool, + const Boundary &) const { Assert (false, ExcPureFunctionCalled()); }; @@ -247,6 +249,7 @@ void FiniteElementBase::fill_fe_values (const typename Triangulation:: template void FiniteElementBase::face_ansatz_points (const typename Triangulation::face_iterator &, + const Boundary &, vector > &) const { Assert (false, ExcPureFunctionCalled()); }; @@ -271,7 +274,8 @@ void FiniteElement<1>::fill_fe_values (const Triangulation<1>::cell_iterator &ce vector > &ansatz_points, const bool compute_ansatz_points, vector > &q_points, - const bool compute_q_points) const { + const bool compute_q_points, + const Boundary<1> &) const { Assert (jacobians.size() == unit_points.size(), ExcWrongFieldDimension(jacobians.size(), unit_points.size())); Assert (q_points.size() == unit_points.size(), @@ -313,6 +317,7 @@ void FiniteElement<1>::fill_fe_values (const Triangulation<1>::cell_iterator &ce void FiniteElement<1>::face_ansatz_points (const typename Triangulation<1>::face_iterator &, + const Boundary<1> &, vector > &) const { // is this function useful in 1D? Assert (false, ExcPureFunctionCalled()); @@ -336,13 +341,15 @@ void FiniteElement<2>::fill_fe_values (const Triangulation<2>::cell_iterator &, vector > &, const bool, vector > &, - const bool) const { + const bool, + const Boundary<2> &) const { Assert (false, ExcPureFunctionCalled()); }; void FiniteElement<2>::face_ansatz_points (const typename Triangulation<2>::face_iterator &, + const Boundary<2> &, vector > &) const { // is this function useful in 1D? Assert (false, ExcPureFunctionCalled()); diff --git a/deal.II/deal.II/source/fe/fe_lib.linear.cc b/deal.II/deal.II/source/fe/fe_lib.linear.cc index 0604e9e7a2..5eafa7b819 100644 --- a/deal.II/deal.II/source/fe/fe_lib.linear.cc +++ b/deal.II/deal.II/source/fe/fe_lib.linear.cc @@ -85,17 +85,19 @@ void FELinear<1>::fill_fe_values (const Triangulation<1>::cell_iterator &cell, vector > &ansatz_points, const bool compute_ansatz_points, vector > &q_points, - const bool compute_q_points) const { + const bool compute_q_points, + const Boundary<1> &boundary) const { // simply pass down FiniteElement<1>::fill_fe_values (cell, unit_points, jacobians, compute_jacobians, ansatz_points, compute_ansatz_points, - q_points, compute_q_points); + q_points, compute_q_points, boundary); }; void FELinear<1>::face_ansatz_points (const Triangulation<1>::face_iterator &, + const Boundary<1> &, vector > &) const { Assert (false, ExcPureFunctionCalled()); }; @@ -246,7 +248,8 @@ void FELinear<2>::fill_fe_values (const Triangulation<2>::cell_iterator &cell, vector > &ansatz_points, const bool compute_ansatz_points, vector > &q_points, - const bool compute_q_points) const { + const bool compute_q_points, + const Boundary<2> &) const { Assert (jacobians.size() == unit_points.size(), ExcWrongFieldDimension(jacobians.size(), unit_points.size())); Assert (q_points.size() == unit_points.size(), @@ -335,6 +338,7 @@ void FELinear<2>::fill_fe_values (const Triangulation<2>::cell_iterator &cell, template void FELinear::face_ansatz_points (const typename Triangulation::face_iterator &face, + const Boundary &, vector > &ansatz_points) const { Assert (ansatz_points.size() == (1<<(dim-1)), typename FiniteElementBase::ExcWrongFieldDimension (ansatz_points.size(), @@ -362,12 +366,13 @@ void FEQuadratic<1>::fill_fe_values (const Triangulation<1>::cell_iterator &cell vector > &ansatz_points, const bool compute_ansatz_points, vector > &q_points, - const bool compute_q_points) const { + const bool compute_q_points, + const Boundary<1> &boundary) const { // simply pass down FiniteElement<1>::fill_fe_values (cell, unit_points, jacobians, compute_jacobians, ansatz_points, compute_ansatz_points, - q_points, compute_q_points); + q_points, compute_q_points, boundary); }; @@ -418,7 +423,8 @@ void FEQuadratic<2>::fill_fe_values (const Triangulation<2>::cell_iterator &, vector > &ansatz_points, const bool, vector > &q_points, - const bool) const { + const bool, + const Boundary<2> &) const { Assert (jacobians.size() == unit_points.size(), ExcWrongFieldDimension(jacobians.size(), unit_points.size())); Assert (q_points.size() == unit_points.size(), @@ -446,12 +452,13 @@ void FECubic<1>::fill_fe_values (const Triangulation<1>::cell_iterator &cell, vector > &ansatz_points, const bool compute_ansatz_points, vector > &q_points, - const bool compute_q_points) const { + const bool compute_q_points, + const Boundary<1> &boundary) const { // simply pass down FiniteElement<1>::fill_fe_values (cell, unit_points, jacobians, compute_jacobians, ansatz_points, compute_ansatz_points, - q_points, compute_q_points); + q_points, compute_q_points, boundary); }; @@ -494,7 +501,8 @@ void FECubic<2>::fill_fe_values (const Triangulation<2>::cell_iterator &, vector > &ansatz_points, const bool, vector > &q_points, - const bool) const { + const bool, + const Boundary<2> &) const { Assert (jacobians.size() == unit_points.size(), ExcWrongFieldDimension(jacobians.size(), unit_points.size())); Assert (q_points.size() == unit_points.size(), diff --git a/deal.II/deal.II/source/numerics/assembler.cc b/deal.II/deal.II/source/numerics/assembler.cc index cf47c72386..616ac3b0ab 100644 --- a/deal.II/deal.II/source/numerics/assembler.cc +++ b/deal.II/deal.II/source/numerics/assembler.cc @@ -51,7 +51,8 @@ AssemblerData::AssemblerData (const DoFHandler &dof, dVector &rhs_vector, const Quadrature &quadrature, const FiniteElement &fe, - const UpdateFields &update_flags) : + const UpdateFields &update_flags, + const Boundary &boundary) : dof(dof), assemble_matrix(assemble_matrix), assemble_rhs(assemble_rhs), @@ -59,7 +60,8 @@ AssemblerData::AssemblerData (const DoFHandler &dof, rhs_vector(rhs_vector), quadrature(quadrature), fe(fe), - update_flags(update_flags) {}; + update_flags(update_flags), + boundary(boundary) {}; @@ -80,7 +82,8 @@ Assembler::Assembler (Triangulation *tria, fe(((AssemblerData*)local_data)->fe), fe_values (((AssemblerData*)local_data)->fe, ((AssemblerData*)local_data)->quadrature, - ((AssemblerData*)local_data)->update_flags) + ((AssemblerData*)local_data)->update_flags), + boundary(((AssemblerData*)local_data)->boundary) { Assert (matrix.m() == dof_handler->n_dofs(), ExcInvalidData()); Assert (matrix.n() == dof_handler->n_dofs(), ExcInvalidData()); @@ -98,7 +101,8 @@ void Assembler::assemble (const Equation &equation) { fe_values.reinit (Triangulation::cell_iterator (tria, present_level, present_index), - fe); + fe, + boundary); const unsigned int n_dofs = dof_handler->get_selected_fe().total_dofs; // clear cell matrix diff --git a/deal.II/deal.II/source/numerics/base.cc b/deal.II/deal.II/source/numerics/base.cc index e5b3cd1875..6e0aff2d48 100644 --- a/deal.II/deal.II/source/numerics/base.cc +++ b/deal.II/deal.II/source/numerics/base.cc @@ -77,7 +77,8 @@ void ProblemBase::assemble (const Equation &equation, const Quadrature &quadrature, const FiniteElement &fe, const UpdateFields &update_flags, - const DirichletBC &dirichlet_bc) { + const DirichletBC &dirichlet_bc, + const Boundary &boundary) { Assert ((tria!=0) && (dof_handler!=0), ExcNoTriaSelected()); system_sparsity.reinit (dof_handler->n_dofs(), @@ -105,7 +106,8 @@ void ProblemBase::assemble (const Equation &equation, right_hand_side, quadrature, fe, - update_flags); + update_flags, + boundary); active_assemble_iterator assembler (tria, tria->begin_active()->level(), tria->begin_active()->index(), @@ -127,7 +129,7 @@ void ProblemBase::assemble (const Equation &equation, // in the docs apply_dirichlet_bc (system_matrix, solution, right_hand_side, - dirichlet_bc, fe); + dirichlet_bc, fe, boundary); }; @@ -158,7 +160,8 @@ void ProblemBase::integrate_difference (const Function &exact_sol dVector &difference, const Quadrature &q, const FiniteElement &fe, - const NormType &norm) const { + const NormType &norm, + const Boundary &boundary) const { Assert ((tria!=0) && (dof_handler!=0), ExcNoTriaSelected()); Assert (fe == dof_handler->get_selected_fe(), ExcInvalidFE()); @@ -187,7 +190,7 @@ void ProblemBase::integrate_difference (const Function &exact_sol { double diff=0; // initialize for this cell - fe_values.reinit (tria_cell, fe); + fe_values.reinit (tria_cell, fe, boundary); switch (norm) { @@ -375,7 +378,8 @@ void ProblemBase::apply_dirichlet_bc (dSMatrix &matrix, dVector &solution, dVector &right_hand_side, const DirichletBC &dirichlet_bc, - const FiniteElement &fe) { + const FiniteElement &fe, + const Boundary &boundary) { Assert ((tria!=0) && (dof_handler!=0), ExcNoTriaSelected()); Assert (dirichlet_bc.find(255) == dirichlet_bc.end(), ExcInvalidBoundaryIndicator()); @@ -387,7 +391,7 @@ void ProblemBase::apply_dirichlet_bc (dSMatrix &matrix, // the lines in 2D being subject to // different bc's), the last value is taken map boundary_values; - make_boundary_value_list (dirichlet_bc, fe, boundary_values); + make_boundary_value_list (dirichlet_bc, fe, boundary, boundary_values); map::const_iterator dof, endd; const unsigned int n_dofs = (unsigned int)matrix.m(); @@ -477,6 +481,7 @@ void ProblemBase::apply_dirichlet_bc (dSMatrix &matrix, void ProblemBase<1>::make_boundary_value_list (const DirichletBC &, const FiniteElement<1> &, + const Boundary<1> &, map &) const { Assert ((tria!=0) && (dof_handler!=0), ExcNoTriaSelected()); Assert (false, ExcNotImplemented()); @@ -487,8 +492,9 @@ ProblemBase<1>::make_boundary_value_list (const DirichletBC &, template void -ProblemBase::make_boundary_value_list (const DirichletBC &dirichlet_bc, +ProblemBase::make_boundary_value_list (const DirichletBC &dirichlet_bc, const FiniteElement &fe, + const Boundary &boundary, map &boundary_values) const { Assert ((tria!=0) && (dof_handler!=0), ExcNoTriaSelected()); @@ -521,7 +527,7 @@ ProblemBase::make_boundary_value_list (const DirichletBC &dirichlet_bc, face_dofs.erase (face_dofs.begin(), face_dofs.end()); dof_values.erase (dof_values.begin(), dof_values.end()); face->get_dof_indices (face_dofs); - fe.face_ansatz_points (tface, dof_locations); + fe.face_ansatz_points (tface, boundary, dof_locations); (*function_ptr).second->value_list (dof_locations, dof_values); // enter into list